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1*  ABSTRACT  .  ' 

X  The  approximation  of  an  analytic  time  function  in  the  leas 
squares  sense  by  sums  of  exponentials  is  considered  from  several 
different  points  of  view.  In  particular,  we  consider  the  determina-- 
tion  of  the  2 n  complex  parameters  s^J  of  the  function  f  (t)  = 

)  exp  (s^t)  so  that  for  a  given  n  and  f(t),  the  value  of  the 

functional  ‘  i  "  -  1  ' 


J  -  J  Cf(t)  -  fa( t) 3  dt 

o  a 

is  minimized.  Equivalently,  the  2n  real  parameters  fa^jb^},  of  the 
Laplace  transfofbrT^* — 

..  .  „  „  .  ,  _  _n-l 


Fa(s)  = 


+  a2  s  +  . . .  +  an  s 
bj  +  bs  s  +  ...  +  bn  s"’1  +  sn 


of  f_(t),  may  be  determined  to  achieve  the  same  minimum  value  of 
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these  two  linear  iterative  schemes  using  several  numerical 
examples . 
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cussed  in  . detail.  Examples  show  the  numerical  difficulties 
due  to  roundoff  errors  that  arise  even  with  the  straight¬ 
forward  methods  available  to  find  the  {o^}.  A  simple 
criterion  for  estimating  these  errors  before  finding  the 
{a^}  is  developed  to  permit  one  to  evaluate  the  feasibility 
of  obtaining  accurate  results  in  any  given  situation. 
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LEAST-SQUARES  approximation  of  functions  by  exponentials 

ABSTRACT 

The  approximation  of  an  analytic  time  function  in  the  leaat-equarea 
sense  by  suns  of  exponentials  is  considered  from  several  different  points 
of  view.  In  particular,  ve  consider  the  determination  of  the  2n  complex 
parameters  of  the  function  fft(t)  ■  exp  (a^t)  *o  that  for 

a  given  n  and  f(t),  the  value  of  the  functional 

00  ? 

J  -  J  lf(t)  -  fa(t)l  dt 

o 

is  minimised.  Equivalently,  the  2n  real  parameters  {a^,b^},  of  the 
Laplace  transfora 

F 

of  fa(t),  may  be  determined  to  achieve  the  same  minimum  value  of  error,  J. 

McDonough's  method  for  finding  the  {s^}  is  derived  by  three  different 
approaches.  In  the  process,  a  new  method  is  developed  vhich  offers  the 
advantages  of  the  earlier  results  achieved  by  McDonough  and  by  McBride, 
Schaefgen,  and  Steiglits.  This  nev  method  reveals  the  link  between  these 
earlier  methods  and  provides  a  standard  for  comparing  these  two  linear 
iterative  schemes  using  several  msMrical  examples. 

The  linear  least-squares  approximation  procedure  in  which  both  n  and 
the  {sk>  (or  {bk})  are  fixed  is  also  discussed  in  detail.  Examples  shew 
the  numerical  difficulties  due  to  roundoff  errors  that  arise  even  with  the 
straightforward  methods  available  to  find  the  A  simple  criterion 

for  estiuating  these  errors  before  finding  the  (e^)  is  developed  to  permit 
one  to  evaluate  the  feasibility  of  obtaining  accurate  results  in  any  given 
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LEAST-SQUARES  APPROXIMATION  OF  FUNCTIONS  BY  EXPONENTIALS 


I.  INTRODUCTION 

1.1  Exponent ial  Representations 

In  approximating  a  function  of  time,  such  at  might  arise  in  control 
or  communication  problems,  one  commonly  uses  a  linear  casibinatlon  of  a  fin¬ 
ite  set  of  simpler  functions.  Exponential  functions  s re  particularly  appro¬ 
priate  for  this  purpose  because  they  have  simple  mathematical  properties. 

It  has  been  demonstrated  in  [l]-[1»]t  that  exponentials  have  very  good 
approximation  properties  for  a  rather  broad  range  of  signal  wave  shapes. 
Furthermore,  in  linear  time-invarient  systems,  the  class  of  exponential 
functions  provides  a  natural  representation  since  the  natural  response  of 
these  systems  is  composed  of  exponential  components.  Another  feature  of 
an  exponential  representation  is  that  there  are  arbitrarily  many  different 
discrete  sets  {expU^t}}  that  are  ccmplet*  over  the  semi-infinite  interval 
with  respect  to  the  L^  norm  (i.e.  the  mean-square  error  in  approximating 
any  piecewise  continuous  f(t)  that  is  square  integrable  over  0  <  t  <  -  by 
the  font  ^  o^exp(sfct)  can  be  made  arbitrarily  small).  This  completeness 
property  is  established  by  Szasx's  form  of  Hunts's  theorem  [3],  which  when 
applied  to  this  exponential  basis  may  be  stated  as  follows:  The  basis 
(expCs^t)}  is  fundamental  with  respect  to  the  L^  norm  over  the  Bead-infinite 
interval  if  and  only  if 


I  - 


Re(sk) 


k-1  !♦ | *k*  ||2 


(1.1) 


Whole  nutbers  in  brackets  refer  to  references  listed  beginning  on  page  81. 


However,  for  practical  work  we  are  not  interested  in  letting  k  approach 
infinity.  Instead,  we  seek  efficient  representation  in  which  k  is  small. 

Of  course,  any  finite  representation  will  necessarily  he  approximate  and 
incomplete.  We  are  interested  in  finding  the  basis  of  lowest  possible 
dimension  thrt  will  lead  to  an  approximation  of  acceptable  accuracy. 

Efficient  representation  will  enable  ua  to  extract  the  information-bearing 
attributes  of  the  signal  with  a  minimum  of  processing.  When  the  interval 
of  approximation  is  finite,  one  can  resort  to  the  discrete  Fourier  series 
since  sines  and  cosines  belong  to  the  class  of  the  exponential  functions. 

But,  despite  the  popularity  of  Fourier  representation,  one  can  often  do 
better  than  this  for  pulse-like  signals  by  using  more  general  exponential 
components.  For  this  reason  exponential  functions  play  an  important  role 
in  s:.gnt,l  representation. 

To  best  approximate  a  signal  by  a  set  of  n  exponentials,  one  must 
determine  the  n  "optimum"  exponents  s^  as  well  as  the  n  amplitudes  o^. 

These  exponents  and  amplitudes  may  be  chosen  to  minimise  the  error  with 
respect  to  some  norm.  Two  popular  norms  are  the  Integrated  squared  error 
(L2  norm) 

J  ■ 

o  k-1  “  “  o 

and  the  Chebyshev  norm  (uniform  norm),  je(t)|,  Hie  former,  often  re¬ 
ferred  to  as  the  least-squares  (or  minimum-error  energy)  criterion  has  been 
studied  extensively  because  it  is  the  most  tractable  mathematically.  For 
given  (s^K  it  i0  easy  in  principle  to  choose  the  for  the  least-squares 
criterion,  since  f&(t)  is  a  linear  function  of  the  {o^}.  However,  practical 
computational  difficulties  exist  because  the  exponential  functions  are  highly 


/  [f(t)  -  l  cl  exp(s.  t)J2dt  *  /  e*(t)  dt 


(1.2) 
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correlated.  As  a  consequence,  solutions  of  the  {a^}  may  be  subject  to 
large  errors  due  to  roundoff  in  the  numerical  computation.  This  difficulty 
will  be  examined  further  in  chapter  II. 

Difficulty  of  a  more  serious  nature  arises  in  finding  the  exponents 
{s^}  for  a  given  f(t)  that  satisfy  the  minimum  error  energy  criterion. 

Until  recently,  only  gradient  methods  were  available,  and  these  frequently 
proved  to  be  quite  unvieldly  for  large  n.  Then  in  1966  McBride,  Steig- 
litz,  and  Schaefgen  [6]  and  in  1968  McDonough  and  Huggi”*  [7]  developed 
two  different  linear  iterative  schemes  which  have  been  found  to  be  quite 
successful  for  determining  the  {s^>  even  for  large  n.  Tvo  natural  ques¬ 
tions  about  these  methods  are  the  following.  First,  how  are  these  meth¬ 
ods  related?  Second,  when  is  it  advent ageous  to  use  one  method  rather 
than  the  other?  This  thesis  provides  answers  to  these  questions  by 
developing  a  new  linear  iterative  method  under  the  least-squares  cri¬ 
terion. 

The  Chebyshev  or  uniform-norm  criterion  has  been  studied  less  than 
the  least-squares  criterion  because  it  is  analytically  more  difficult. 
Apparently,  not  much  has  been  done  with  thiB  criterion  to  date  but, 

Tang  [8]  has  shown  how  the  {a^}  may  bfe  determined  provided  all  the  {s^} 
are  real  and  distinct.  So  far,  it  appears  that  the  only  way  to  find  the 
exponents  {s^}  for  the  Chebyshev  criterion  is  by  slowly  converging  grad¬ 
ient  methods. 

1.2  Some  Known  Methods  of  Approximation  by  Exponentials 

1.2.1  Non-Qptimal  Approximation  -  Prony's  Method  and  Pade  Approximants 

Two  simple,  but  often  successful  ways  of  obtaining  an  approximation 


to  a  function  by  sums  of  exponentials  use  Prony's  method  and  Pade  ap- 
proximants.  Neither  results  in  an  optimal  approximation  with  respect 
to  the  Lg,  uniform,  or  any  other  norm,  but  they  do  provide  tvo  quick 
and  straightforward  ways  of  obtaining  approximations  that  are  usually 

"fairly  good".  In  Pade  approximants  one  matches  the  rational  function 

i;  ,  k-i 

F  (a)  «•  b  .  *  1  (1.3) 

a  rn+1  .  n-i  D\s)  n+1 

ik-1  “k  3 

to  the  desired  function  F(s)  (the  Laplace  transform  of  f(t)  )  by  adjust¬ 
ing  the  {a^.b^}  8uch  that  Fft( a )  will  have  the  same  power  series  as  the 
power  series  expansion  of  F(s)  through  the  m+n  power  where  m<n.  That 
is,  the  seminorm 

| 1 F(s )— F  (a ) | |=|F(0)-F  (0)|+|F'(0)-F'(0)|+...+jFm+n(0)-F  “+n(0)|  (l.U) 

A  A  ft  ft 

is  made  zero.  The  real  merit  of  the  Pade  method  is  the  computational 
ease  with  which  the  may  be  found.  Finding  the  (b^)  involves 

solving  n  linear  equations  in  n  unknowns.  Once  the  (b^)  are  determined, 
the  {a^}  are  similarly  found  by  evaluating  another  linear  system  of  m 
equations.  These  are  explicit  equations,  not  simultaneous  for  the  (a^). 

To  write  Fa(s)  as  a  sum  of  exponentials,  a  partial-fraction  expan¬ 
sion  must  be  performed  which  requires  finding  the  roots  of  the  nth 
degree  polynomial  D(s).  Kautz  [9]  and  Mathers  [10]  have  used  the  method 
in  designing  circuits  to  approximate  specified  transient  responses. 
Teasdale  [ll]  first  applies  the  transformation  z»(l-s)/(l+s)  to  ob¬ 
tain  an  "indirect  Pade  approximant"  matching  a  power  series  in  i  instead 
of  s.  (Actually,  since  z*0  implies  s*l,  this  is  matching  the  power 


series  stout  the  point  s»l  instead  of  s»0.)  The  procedure  developed 
will  be  different  from  the  direct  Pade  approximate t  with  generally  small¬ 
er  error  but  at  the  expense  of  more  computation. 

Another  simple  way  of  approximating  a  function  by  sums  of  ex¬ 
ponentials  is  a  method  first  used  by  Prony  in  1795.  This  procedure 
was  first  applied  to  network  synthesis  problems  by  Tuttle,  Carr,  and 
Kautz.  A  detailed  discussion  of  the  method  and  its  refinements  is 
given  in  McDonough's  thesis  [ 1 ] .  The  principle  behind  the  method 
originates  from  the  fact  that  if  a  waveform  is  indeed  composed  of 
exponentials,  viz. 


n 

f(t)  ■  l  cl  exp  (s^t)  Re{8k>  <  0 

k»l 


then  f(t)  will  be  the  solution  to  some  homogeneous  differential 
equation  of  the  nth  order, 


n  ..i- 

z 

i*o  1  dt 


B 


(1.5) 


(1.6) 


Provided  one  could  find  the  coefficients  (B^)  of  this  equation,  the 
{ 8^}  could  then  be  obtained  by  evaluating  the  n  roots  of  the  poly- 
nominal  BjS^O.  Our  task  then  is  to  find  the  B^  appropriate  to 
a  given  f(t).  Then,  the  (sk>  which  satisfy  the  differential  equation 
may  be  used  to  construct  an  exponential  basis  for  f(t).  If  the  signal 
is  noisy  or  is  not  exactly  the  sub  of  n  exponentials,  the  left  hand 
aide  of  equation  (1.6)  cannot  be  made  zero  regardless  of  the  choice 
of  (B. }  and  there  will  be  a  residual  LJJ  B,  (d^f/dt*)  ■  e  (t). 

X  1"0  i  p 


-11- 


Since  Bq*1,  equation  (1.6)  nay  be  written  as 

Su>  ■ f(t)  *  X  B*  iftil  (1.T) 

dt 

Hien,  one  simply  chooses  the  set  of  (B. )  to  minimize  this  e  (t)  in  the 

1  P 

least-squares  sense,  thus 

E  •  /  [c  (t)]2  dt.  (1.8) 

o  y 

Minimizing  E  over  the  coefficients  re8U^8  n  linear 

simultaneous  equations 


/  f*i}f  dt  +  l  B  /  f^  f(k)  dt  -  0 
o  k«l  o 


i  *  1,2,.. .  ,n. 


(1.9) 


However,  the  matrix  elements 

gik  *  7  fU)  fU)  dt  (1.10) 

o 

will  not  exist  unless  f  is  of  at  least  class  c“.  If  the  differential 
equation  is  first  integrated  n  times  the  corresponding  new  elements 
will  exist  for  any  piecewise  continuous  function  with  finite  energy 
but,  this  initial  integration  should  be  performed  the  least  mmiber  of 
times  to  assure  the  existence  of  (1.10)  since  it  has  a  tendency  to 
destroy  signal  information.  Fortunately,  the  matrix  elements  have 
certain  recursion  relations  which  make  it  necessary  to  calculate  the 
gkk  on^*  irony’s  method  yields  only  the  frequencies  (s^)  but,  the 
smplltude  coefficients  {o^}  may  subsequently  be  found  with  little 
difficulty  (as  discussed  in  the  next  chapter).  It  should  be  emphasized 
that  Prony's  method  does  not  lead  to  the  optimia  least-square  approx- 
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imation  (unless  f(t)  is  exactly  the  sun  of  n  exponentials)  since 

e  (t)  is  not  identical  to  f(t)-f  (t). 

P  * 

1.2.2  Optimal  Approximation  in  the  Least-Squares  Sense  by  Exponentials 
The  conditions  for  optimal  exponential  approximation  of  a  func¬ 
tion  f(t)  with  respect  to  the  L£  norm  over  the  semi-infinite  interval 
are  described  compactly  oy  the  equations  of  Aigrain  and  Williams  [12]. 
Although  theoretically  attractive,  these  nonlinear  transendental 
equations  are  computationally  undesirable  and  algebraic  solution  is 
seldom  possible  even  when  the  Laplace  transform  of  f(t)  is  known  in 
closed  form.  Most  often,  these  equations  are  solved  by  a  gradient 
or  some  other  iterative  method. 

In  chapter  III,  it  will  be  shown  that  by  suppressing  the  amplitude 
coefficients  (a^)  one  may  write  the  Integrated  square  error,  J,  as 

m 

J  -  /  f*(t)  dt  -  F  H’1  F  (1.11) 

o 

where  is  the  inverse  of  the  generalised  Hilbert  matrix  and  £  is  a 
column  matrix  (^[F(s*)  ,F(a") ,. . .  ,F(s*)  ] ).  Equation  (l.ll)  is  a  com¬ 
pact  mathematical  expression  for  the  mean-square  error  but,  for  large 
n,  (say  n>£)  it  is  very  sensitive  to  the  variation  of  the  (s^J  and 
finding  the  minimum  J  by  the  usual  gradient  methods  may  be  ineffec¬ 
tive.  The  two  mu.h  more  effective  ways  of  solving  the  Aigraln-Williams 
equations  for  large  n,  which  have  been  recently  developed,  shall  new 
be  briefly  discussed. 

The  method  of  McBride,  Scbaefgen,  and  Steiglits,  the  first  of  the 
linear  iterative  methods  mentioned  earlier,  introduces  an  approximate 
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error 


D,(a)  N  (a) 

E.(>)  ■  «•>  -  s^Itt 


with 


r  (.)  -SW 

a'  '  D(s) 


w 


...+a> 

n 


n-1 


.  . .  , ,  n-1,  n 

b-+b0B*...+b  a  +s 
id  n 


(1.12) 


(1.13) 


where  the  subscript  J  refera  to  the  iteration  number.  The  previously 

th 

computed  coefficients  of  are  regarded  as  fixed  during  the  J 
iteration.  By  this  simple  tactic,  the  error  is  linearized  in  terns  of 
the  unknown  coefficients  { a^,b^}  of  the  numerator  polynomials  and 
D..  The  primary  difficulty  with  this  method  is  that  the  approximate 

V 

rather  than  the  true  error  is  be  inf;  minimized.  Hence,  the  iterative 
scheme  does  not  converge  to  the  true  optimum.  To  circumvent  this  diffi¬ 
culty,  McBride  et.  al.  introduce  a  different  "Mode-2  Iteration"  which 
does  converge  to  the  true  minimus  but  more  slowly  than  one  would  hope. 

The  requirement  for  using  two  different  iteration  schemes  also  adds  extra 
complexity  to  the  McBride  method. 

The  difficulties  of  the  McBride  method  are  avoided  in  the  linear 
Iterative  scheme  devised  by  McDonough  and  Huggins.  Here,  the  2n 
Aigrain-Williasm  equations  are  first  reduced  to  a  set  of  n  equations 
involving  the  {s^>  only.  This  was  done  by  regarding  F( a )  as  a  signal 
in  a  vector  apace  and  showing  that  a  necessary  condition  for  the  (s^) 
to  be  optimum  is  that  F(s)  be  orthogonal  to  the  space  spanned  by  +^(s) 

1  -  1,2,...  ,n  where 


i-1 


♦*(•) 


H(s 


( ses  J 


''-Vi 


¥  S-B 


i  k-1 


(1.1*0 
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with 

This  orthogonality  condition  may  be  written  as 

!  Fl->  5^7  *  0  1  ■  1.2 . »• 


-J* 


(1.15) 


(1.16) 


The  linear  iterative  scheme  described  by  McDonough  is  obtained  by  re¬ 
placing  H(s)  with  the  new  operator 
n+1  .  . 

Ha(s)  -  l  b1(-si)l*i  /  D(s)  ,  bn+1  -  1  (1.17) 

The  resulting  iterative  method  is  similar  to  the  one  described  by  McBride. 
All  these  optimum  least-squares  methods  will  be  discussed  more  fully  in 
chapter  Ill. 

1.3  Brief  Discussion  of  Previous  Methods 

Prony's  method  and  the  method  of  Pade  approximsnts  have  two  things 
in  cannon.  First,  each  requires  the  solving  of  a  system  of  linear 
simultaneous  equations.  Second,  to  find  the  approximate  { s^} ,  one  must 
evaluate  the  roots  of  an  nth  degree  polynomial.  Each  method  uses  the 
application  of  these  two  operations  only  once.  Hence,  each  is  useful 
in  that  it  provides  a  rapid  way  of  obtaining  an  approximation  to  a 
desired  vaveform.  To  improve  these  initial  approximations  or  to  make 
them  optimal,  either  of  the  linear  iterative  schemes  may  be  used. 

These  linear  iterative  methods  also  require  solving  a  system  of  linear 
simultaneous  equations  and  finding  the  roots  of  an  n*h  degree  polynomial 
for  each  iteration.  It  vill  be  shewn  in  the  thesis  that  the  method  of 
McDonough  is  the  better  way  of  improving  the  approximation. 


Regardless  of  the  initial  starting  point  in  the  approximation, 
Mode— 1,  Mode-2,  and  McDonough's  method  all  converge  to  the  optimal 
solution  in  one  step  if  the  function  f(t)  is  exactly  coir.poaeil  of  n 
exponentials.  This  suggests  that  any  of  the  linear  methods  will  con¬ 
verge  rapidly  to  the  optimum  when  the  signal  is  "nearly  exponential". 
I.l»  Plan  of  the  Thesis 

It  is  veil  known  hew  to  find  the  amplitude  coefficients  for  a 
least-square  approximation  to  a  function  on  a  fixed  or  known  basis. 
However,  computational  difficulties  arise  when  the  basis  elements  are 
highly  correlated.  In  chapter  II  a  closed-form  expression  is  developed 
for  the  inverse  of  some  Gram  matrices  that  occur  in  least-square 
theory.  These  new  expressions  can  help  to  reduce  roundoff  errors 
in  computing  the  amplitude  coefficients.  In  particular,  the  exponen¬ 
tial  basis  is  studied.  An  explicit  inverse  for  the  generalized  Hil¬ 
bert  matrix,  the  Gram  matrix  for  an  exponential  basis,  was  published 
in  a  French  Journal  in  i960  [lUJ.  Although  this  result  is  quite  use¬ 
ful  in  least-square  representation  by  exponentials  and  polynomials, 
it  appears  to  have  remained  unknown  to  the  English  literature.  Its 
use  is  fully  explored  in  regard  to  finding  the  amplitudes  as  veil  as 
the  {s^}. 

The  successful  methods  of  McDonough  and  Huggins  and  of  McBride, 
Schaefgen,  and  Stelglitz  are  discussed  in  detail  in  chapter  III.  A 
nev  scheme  is  developed  which  for  the  first  time  enables  one  to  make 
a  meaningful  comparison  of  the  methods.  In  chapter  IV,  several 
numerical  examples  compare  the  convergence  properties  of  the  linear 
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Iterative  methods.  Finally,  in  chapter  V,  the  new  method  is  extended 
to  deal  vith  imperfectly  known  signals  or  sampled  data. 
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II.  DETERMINATION  OF  THE  AMPLITUDE  COEFFICIENTS  IN  LEAST-SQUARE  APPROX¬ 
IMATION  OF  FUNCTIONS  BY  EXPONENTIALS  AND  OTHER  COWO.LY  USED  BASIS 
FUNCTIONS 


Suppose  ( t ) ,  x2 ( t ) ,  ....  xn(t)  denote  a  finite  sequence  of  con¬ 
veniently  chosen  functions  defined  over  some  continuous  interval  (a,b) 
of  t.  Let  f^(t)  ■  1^2.  be  911  aPProxiBfttion  t0  the  function 

f(t).  One  problem  is  to  find  the  such  that  /^lf(t)  -  f^(t)|2  dt  ■  J 
is  a  minimum.  The  standard  least-squares  procedure  yields  the  following 
equations: 

aj  n 

~  “  0  °r  l  VS  *  fi  1  "  *»2 . n  (2<1) 

i  f  1  J  3 

where  ■  /“  x^tix^t)  dt  »  g*^  are  the  elements  of  the  Grsm  matrix 
0  end  »  /”  f(t)x*(t)  dt  are  the  elements  of  the  column  F.  Then  the 
best  fitting  amplitude  coefficients  are  given  by  the  column  matrix  A, 
where 

A  •  G-S.  (2.2) 

Theoretically,  this  algorithm  presents  no  difficulty  provided  the  x^  are 
linearly  independent.  If  they  are  not,  the  matrix  G  will  be  singular 
but  the  same  minimum  error  J  can  be  obtained  with  a  set  of  less  than 
n  of  the  x^  that  are  independent.  When  the  x^  are  highly  correlated, 
but  independent  (exponentials  for  example),  the  matrix  G  is  "ill- 
conditioned"  and  costputational  difficulties  arise  in  finding  the 
inverse  accurately  for  any  sixeable  n,  as  evidenced  in  [l9)*[2l]. 

This  difficulty  is  saietimes  reduced  by  Introducing  a  new  basis  of 
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orthcnormal  functions,  which  are  linear  combinations  of  the  original 
i>aals  functions  x^(t)  and  span  the  ssme  space.  However,  these  ortho¬ 
normal  functions  no  longer  poaess  the  simple  properties  of  the  origi¬ 
nal  basis  so  this  approach  is  not  a  cure-all.  A  method  for  finding 
is  still  needed. 


II. 1  Inverse  of  the  Gram  Matrix 

Let  41(t),t2(t),...,#n(t)  be  a  set  of  orthonormal  functions, 
which  may  be  determined  from  the  x^(t)  by  the  Gram-8chmidt  procedure. 
That  is, 

i 

♦ilt>  ’  I  cuV‘>  i  ■  1,2,*..  ,n  (2.3) 

and  c^  cannot  be  zero  if  the  x^  are  linearly  independent.  Written  in 
matrix  notation,  #(t)  ■  CX(t),  where  C  is  a  nonsingular  n  X  n  lever 
triangular  matrix  and 


/  ♦i(t)^J(t)  dt  -  6±J 


(2.M 


Using  Dirac  notation,  let  X|  denote  the  column  of  basis  elements  and 
|  )C  its  adjoint.  Then 


From  (2.3) 


Xj  X  -  G 

(2.5) 

il  •  C  xj 

(2.6) 

II  •  £  xil 

(2.7) 

II-  li  £ 

(2.8) 

lefts  that  1*  the  adjoint  of  £  which 
real,  then  C  is  the  transpose  of  C. 
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x|x  -  c  ♦  c  , 


(2.9) 


(identity  matrix), 


Hence 


G  -  C"1?"1 


(2.10) 


or  the  final  result 
G  «  C  C. 


(2.11) 


Equation  (2,11)  is  a  useful  result  in  two  ways.  Firsts  it  nay  be  used 
to  find  explicit  expressions  for  the  inverses  of  some  Gram  matrices 
that  are  ill-conditioned.  Second,  it  can  sometimes  be  used  to  construct 
an  orthonormal  basis  by  simple  inspection.  The  second  application  is 
not  one  of  the  main  goals  of  the  thesis  and  therefore,  is  discussed 
in  Appendix  A.  The  first  application  will  now  be  demonstrated  by  find¬ 
ing  the  inverse  of  the  Hilbert  matrix. 


II. 2  The  Generalized  Hilbert  Matrix 


The  generalized  Hilbert  matrix  is  the  n  X  n  Hermitian  matrix  with 


elements  h. 


-  (v**rl  or 


i/(Sl+S*)  l/U^+s*)  ...  l/(SjL+s*) 


l/(s£+s*)  l/Ug+Sg)  ...  1/(s2+b*) 


(2.12) 


|_1/(VS1)  1/(V82>  — 

where  the  si  are  n  complex  scalars  and  s^  4  s^  if  i  4  J,  and  si  4  0. 
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The  Hilbert  matrix  la  diacuaaed  extensively  in  the  literature  [21]- 
[2b].  The  inverse  of  this  matrix  is  shewn  to  have  for  its  elements 

<V.*)(y.*)  r  “  (V.*n  f  -A- 

kyi  ktf) 

Pr->of :  Let  x^(t)  ■  exp  (♦s^t) ,  Re(s^  <  0),  then  h^  * 
x^(t)  x*(t)  dt.  As  in  (2.3) «  1st  the  orthonormal  functions  bs  given  by 
♦  (t)  ■  £  X(t),  where 

'  .  ^  . 


cij  * 


ut*v 


(2.1b) 


The  t^(t)  are  known  to  be  orthonoraal  fran  Kautz's  method  [ 13 J »  [18]. 
The  Laplace  transform  of  ^(t)  is 


(-8.-8*)^^  -r-r-  (s4s£) 

'  Ts«'.'«) .  TT  1ZTT  ■ 

1  Lk-1  k  J 


The  c . .  in  (2.1b)  correspond  to  the  residues  of  this  transform. 
From  (2.11) 

-1  ?  « 

hiJ  ckickJ 

max  (i„J)  - 

since  c.  ■  0  if  i  <  m.  From  (2, lb)  and  (2.l6) 
im 


—  1  u 

h .  ■  c*  c 

in  ni  nn 


(si+s*)( 


ffr  (V*)1  ffr  (*k-‘n)l 

F?  Li  TW  '  i-i  TWJ 


(2.15) 


(2.16) 


(2.17) 


Because  of  the  symmetry  in  the  original  matrix,  if  cfln  is  replaced  by  c^ 
in  (2.17),  the  formula  must  hold  for  the  general  term  h"^  (since  the 
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order  of  Bi*82»“*»8n  in  1L  ctn  be  changed  without  affecting  the  form 
of  the  equations)  and  (2.13)  is  proved.  In  the  special  case  that  all 
the  s^  are  real,  the  formula  reduces  to 

f-TT 

yh  (2.18) 


,.1  f-TT  <Vi>1  f-TT  (V', 

ij  will 


kyi 


kyj 


This  result  agrees  with  Oastinel  [lU]  who  found  this  expression  for  the 
inverse  of  a  generalized  Hilbert  matrix  by  a  rather  tedious  application 
of  Lagrange's  interpolation  polynomial.  Appendix  B  gives  another  inter¬ 
esting  explicit  inverse  using  Laguerre  functions. 


II. 3  Roundoff  Errors  in  the  Amplitude  Coefficients  Using  a  Fixed 
Exponential  Basis 

Let  f(t)  be  a  piecewise  continuous  real  function  having  finite 

energy  in  the  semi-infinite  interval,  i.e. 

00 

/  f*(t)  dt  <  -.  (2,19) 

o 

We  wish  to  find  the  amplitude  coefficients  that  will  minimize  the 

mean- square  error, 

"  n  2 

J  ■  /  [f ( t)  -  l  exp(sRt) ]T  dt  (2.20) 

o  k»l 

for  a  specified  set  of  exponential  functions,  having  {s^}  with  negative 
real  parts.  From  (2.1),  the  simultaneous  equations  for  determining 
the  (o^)  are+ 

+  Since  f(t)  is  real,  the  must  occur  in  complex  conjugate  pairs. 
Hence,  there  is  no  loss  in  generality  if  every  s*  is  replaced  by  s^ 
in  (2.21)  and  to  simplify  the  typography  this  will  be  done  throughout 


the  remainder  of  the  thesis. 
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lASj+S*)  l/Uj+B*) 

1/(b2+§*)  1/(b2+b*) 

•  • 

•  • 

•  • 

l/(s  +S*)  l/(s  *S?) 
n  1  n  d 


1/1 V**) 

• 

l/(a  *s*) 
n  n 


where  F(b)  is  the  Laplace  transform  of  f(t).  In  matrix  form  these 
equations  are  F  ■  H  A,  where  H  is  the  generalized  Hilbert  matrix, 
and  their  solution  is  A1  H**1  F.  However,  the  Hilbert  matrix  is 
notoriously  ill-conditioned  and  computation  of  H-^  by  any  of  the 
standard  methods  (Gauss -Jordan,  Seidel's  method,  method  of  Crout, 
etc.)  encounters  serious  roundoff  difficulties  for  n  greater  than 
5  or  so,  even  when  double-precision  arithmetical  operations  are  used. 

The  rapid  growth  of  roundoff  errors  with  increasing  n  may  be 
demonstrated  by  canparing  the  inverse  of  H  (for  s^i  i*l,2,. , .  ,n 
with  n*5  and  7)  calculated  by  the  explicit  formula  (2.13)  with  the 
inverse  obtained  by  the  method  of  Crout  [15 1.  All  calculations  were 
made  by  an  IBM  709^  having  approximately  8  significant  decimal  digits. 
Table  2.1  shows  that  only  3-place  accuracy  is  attained  in  many  of  the 
elements  of  the  inverse  matrix  for  n*5t  and  a  complete  loss  of  sig¬ 
nificant  results  for  n*7  using  Crout 's  method.  For  still  larger  n 
the  results  are  meaningless.  On  the  other  hand,  the  explicit  formula 
achieves  7-place  accuracy  (for  both  n"5  and  7).  (This  was  validated 
by  double-precision  calculations.)  Since  a  detailed  analysis  of 
roundoff  errors  arising  in  inversion  of  matrices  on  computers  is 
given  in  references  [l6]  and  [17],  this  topic  will  not  be  discussed 


further  here.  Moreover,  the  explicit  formula  for  inverting  the  Hilbert 
matrix  has  so  reduced  these  errors  in  finding  the  inverse  as  to  make 
them  insignificant  for  modest  n. 


Table  2.1 


Inverse  of  the  Hilbert  Matrix  h. ■l/(i+j)  for  n«5  and  7  by  Method 
of  Crout  and  by  Explicit  FormulaJon  Computer  Having  8  Significant 
Figure  Accuracy. 
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Explicit  Formula 

n«5  (Exact  Inverse) 
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Method 

of  Crout  n*5  (3  Significant 

Places ) 
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Explicit  Formula 

n*7  (7  Significant 
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Method  of  Crout  n«7  (No  Significance) 
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Another  source  of  roundoff  error  occurs  in  the  cosputation  of  the 

{o^}  when  H*1  is  multiplied  by  £.  From  (2.21)  and  the  explicit  formula 

•1  th 

for  H  ,  the  k  amplitude  coefficient  can  be  expressed  as 


-2k- 


a  ■  -  1*8.  T? 
k  k  k 


n  s  T? 

i  r+rK-V 

i-1  3i  *k  lJ 


where 


T<V%> 


■TTr8-* 

1 


7 


(2.22) 


(2.23) 


m*l 

B^k 


The  estimation  of  a  roundoff  error  in  evaluating  equation  (2.22)  may 
be  illustrated  by  considering  the  approximation  of  a  square  pulse, 
f(t)  ■  u_j(t)  -  u^ft-l),  where  u^(t)  is  the  unit  step, 

u^t)  ■  1  t  >  C 

■  0  t  <  0 

and  thus 

F(s)  -  (l-e"8)/s.  (2.21») 

Again  let  s^«i  i*l,2,...,n.  Columns  1,  2,  and  6  of  Table  2.2  sumnarixe 
the  results  using  equation  (2.22)  with  n»5,7»and  9.  (The  error  is  es¬ 
timated  by  comparing  these  results  with  those  obtained  by  double- 
precision  calculations.)  Notice  that  and  t”  exhibit  nearly  the 
same  order  of  magnitude  for  almost  all  n  and  k.  This  implies  that 
the  sun  of  the  n  terms  within  the  brackets  of  (2.22)  must  be  roughly 
of  unit  magnitude.  Also,  all  of  the  Ff-s^)  are  less  than  1.  (In 
chapter  V  it  is  shown  that  for  any  normalised  f(t),  | r(— s^) 
(2Re{-sk)]"1^2. )  Therefore,  the  ntaber  of  significant  decimal 
digits  lost  in  each  of  the  e^  computed  by  this  method,  which  forms 
small  differences  of  very  large  maters,  may  be  expected  to  be  the 
same  as  the  number  of  places  to  the  left  of  the  decimal  point  in 
the  largest  t£.  In  the  example  provided  by  Table  2.2,  for  n«5 , 
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has  the  largest  magnitude  of  315 .  Thus,  a  loss  of  about  three  sig¬ 
nificant  places  may  be  expected;  for  n»9t  T^  ■  -210210  indicating  a 
loss  of  six  significant  places  in  each  o^.  These  predictions  agree 
with  the  actual  accuracies  obtained  in  Table  2.2, 

Tsble  2.2 


Amplitude  Coefficients  of  exp(-kt)  k«l,2,...,n  for  a  Least-Square 
Approximation  of  the  Square  Pulse 


cl  From 

4.(2.221  irr°r 

Eq.f2.32)  Error 

a,  Double- 
Precision 

* 
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n»5  Loss  of  3-Places 
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-120. 6H2 
1012.112 
-  16  7.217 

-0.006 
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n«7  Loss  of  5-Places 
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?* 109,999 

n-9  Loss  of  7-Places 


A  procedure  that  is  often  used  to  avoid  the  inversion  of  the 
ill-conditioned  aatrix  of  equation  (2.21)  is  achieved  by  introducing 
orthogonaiixed  combinations  of  the  original  exponentials.  The  orthog- 
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onalizing  procedure  may  be  implemented  in  practical  filters  by  a 

method  due  to  Kautz  [13],  [18]  which  it  based  on  the  traditional 

Gram-Schmidt  procedure  applied  to  exponential  functions.  Rewriting 

(2.1k) ,  the  orthonormal  functions  are 
i 

♦jU)  ■  l  cik  exp(skt)  (2.25) 


where 


ik 


(-1) 


i+1 


(ak^a*)/-(Sj-»s^) 


V'k 


i  >  k 


(2.26) 


Then 


n 

f(t)  -  l  d.  ♦v(t).  (2.27) 

m  k«l  *  “ 

As  is  well-known  from  the  theoxy  of  orthonomal  basis,  the  expansion 
coefficients 

dj  ■  /  f(t)  #k(t)  dt  (2.28) 

o 

automatically  yield  a  least-squares  fit.  Equation  (2.28)  may  be  ex¬ 
pressed  in  matrix  form  as 

(2.29) 

or  D  ■  C  F.  When  the  signal  coordinates  on  the  orthonormal  basis  are 
-*  as  “ 

transformed  to  find  the  coordinates  on  the  original  exponential  basis 


•—  — 

dl 

c^  0  0...  0 

F(-s1) 

d2 

• 

a 

— 

C21  C22  0  0 

•  •  • 

•  •  e 

• 

• 

.  «« 

e  •  • 

C  .  C  -  C  *s  e  e  C 

nl  n2  n3  nn 

• 

one  gets 


°1  "  cudl  *  c21dl  +  C31dl  +  ••• 


C22d2  *  C32d3  +  *•' 


(2.30) 


°33d3  +  *** 


Equation  (2.30)  is  easily  verified  using  (2.2),  (2.1l)  and  (2.29). 

Hence ,  explicit  equations  for  the  (o^)  can  be  obtained  in  tvo 
simple  steps  by  combining'*'  (2.29)  and  (2.30), 

A  -  C  (C  F)  (2.31) 

or 


I  cik  l  cij  F(‘8J) 

i«k  |_J-1  J  J  „ 


(2.32) 


To  minimize  the  nisnber  of  arithmetical  operations  required  in  evaluat- 
ing  A  *  £  g  R,  the  product  £  F  should  be  formed  first.  This  requires 
n(n+l)  multiplications,  whereas  if  C  C  is  formed  first,  roughly 
n  /2  multiplications  are  needed  to  find  all  the  {a^}. 

Although  A  may  be  evaluated  by  either  equation  (2.22)  or  (2.32) 
(which  are  theoretically  equivalent),  equation  (2.22)  is  computa¬ 
tionally  preferred  for  three  reasons.  First,  it  is  much  simpler, 
requiring  only  a  tingle  simulation  rather  than  the  double  simulation 


|  ^  .  . 

Of  course,  if  the  £  £  in  (2.31)  is  combined  and  simplified,  cme 

obtains  the  explicit  expression  for  the  inverse  of  the  generalized 
Hilbert  matrix  obtained  earlier.  Apparently,  this  way  of  finding 
the  inverse  of  a  Gram  matrix  has  not  appeared  previously  in  the 


literature 
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of  (2.32).  Second,  only  n  are  needed  in  the  first  method,  whereas 
in  the  second,  n(n*l)/2  different  c^j  must  be  evaluated.'*’  Third, 
the  method  of  equation  (2.22)  provides  a  simple  estimate  of  the 
number  of  significant  places  that  will  be  lost  due  to  roundoff  even 
before  the  actual  computation  is  made. 

By  observing  the  magnitudes  of  the  l£,  in  Table  2.2,  we  have 
already  noted  that  the  percent  roundoff  error  corresponds  to  the 
magnitude  of  the  largest  T^.  Table  2.2  shows  that  the  accuracy  of 
either  method  is  about  the  same,  so  the  choice  rests  entirely  on 
which  offers  the  greatest  computational  advantage:  this  is  the 
method  of  equation  (2.22). 

Thus  far,  we  have  presented  two  methods  for  determining  the 
amplitude  coefficients.  For  single-precision  computation  size¬ 
able  numerical  errors  arise  in  both  methods  for  n  greater  than  U, 


t 

Some  simplification  is  possible  in  evaluating  these  c^  by  making 
use  of  a  recursion  relation  which  requires  the  calculation  of  only 


the  n  c^,  all  other  quantities  being  obtained  from  these.  The  re¬ 


lation  obtained  from  (2.23)  and  (2.26)  is 


'i+1 


V’k 


i**T 


vr*k 


'ik 


l  >  k 


Even  using  this  recursion,  the  computation  of  the  c^  alone  requires 

at  least  as  much  work  as  the  tT. 

k 
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and  for  n  greater  than  15  (maximum  t£  >  10"°  with  s^»i)  even  double¬ 
precision  arithmetic  may  not  be  adequate.  There  remains  a  real  need 
for  further  detailed  study  of  the  computational  aspects  of  these 
methods. 


11.4  The  Vandermonde  Matrix 

The  Vandermonde  matrix  arises  in  many  branches  of  applied  mathe¬ 
matics.  In  control  theory  one  encounters  the  equation  t  (t)  » 

A&(t)  ♦  D  m(t)  [27]  •  which  may  be  simplified  by  transforming  +he 
state  vector  X  to  Y  ■  V’1  X  where  V  is  the  Vandermonde  matrix. 

In  nunerical  interpolation  by  polynomials  of  a  function  defined 
by  a  set  of  n  ordered  pairs  of  real  or  complex  mmbers  (s^jZ^) 
with  all  the  s^  distinct,  one  seeks  a  unique  polynomial 

N(s)  ■  a^  ♦  ag  s  ♦  ...  ♦  aQ  s0”1  (2.33) 

for  which 

N(sk)  *  ^  k«  1,2,... ,n.  (2.34) 


The  conditions  (2.34)  form  a  system  of  n  linear  equations  in  the 
Sj,  coefficients  of  the  polynomial. 


”  2  n-l~ 

"  — ■ 

m 

Is.  S.  ...  8 

a. 

11  1 

1 

1 

.  2  n-1 

1  *2  *2  •  •  • 

*2 

*2 

M 

...  • 

e 

e 

e 

• 

.  •  •  •  • 

• 

• 

,  2  n-1 

lS  S  .  .  .  6_ 

a 

a 

n  n  m 

_  n„ 

(2.35) 


Hie  matrix  of  this  system  is  named  after  Vandermonde  and  is  shewn  to 
be  non-singular  provided  all  the  s^  are  distinct  [28]. 


30- 


The  Vandermonde  matrix  also  arises  in  least-square  approximation 


using  exponentials  over  the  semi-infinite  interval  as  we  now  show. 

By  solving  equation  (2.21),  one  obtains  the  best  fitting  approximation 


F  (a)  -  — ±-+  ...  ♦  -2- 

a  s-s.  s— s_  s-s 

i  t  n 


(2.36) 


which  has  the  properties  ennunerated  by  Aigrain  and  Williams  [12],  that 


F(-s  )  ■  F  (-s.)  k  ■  1,2, ...,n. 


* '  a'  k/  i*  •  •  I** * 

Equation  (2.36)  may  also  be  written  as  the  rational  fraction 

a,+a_s+...+a  a"”1 
F  (s)  -  -i-2 - S - - - 

&  b,+b0s+...+b  8Xl"’1+sn 

i  c  n 


(2.37) 


(2.38) 


[fl-s.)(s-8,J...Ts-s  )  D(s)  * 
1  e.  n 


(2.39) 


When  the  {a^}  are  optimally  chosen,  equation  (2.37)  is  satisfied.  Then, 
D(-sk)  F(-sR)  «  F^-s^)  *  Nl“8k)  k«l,2,...,n<  (2.1*0) 

In  matrix  fora, 

(ol-).)  F(-.J|  fl  (-.j) 


D(-.2)  r(-.2)  .  1  (-2>  <-.2>:..l-»2>  *2 


(2.1*1) 


Lc(-n»  F<-v|  L1  <"»)  (-v-~  ‘-.HLsJ 

Equation  (2.1*1)  also  exhibits  the  Vandensonde  matrix  V  with  elements 


V  m 

ij 


(— s^)1^  i , J*l,2 , •  • .  ,n . 


(2.1*2) 


To  solve  (2.1*1)  for  the  coefficients  of  the  numerator  polynomial, 
a  closed-form  expression  for  the  inverse  of  the  Vandermonde  matrix 


is  needed.  An  explicit  expression  for  the  Inverse  of  this  important 
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matrix  is  given  by  Tou  [29]  and  the  result  may  be  summarized  by  the 
following  theorem: 

Theorem  1^  -  Let  V  be  the  Vandermonde  matrix  with  elements  v^*(-s^)^"* 
and  V  1  be  its  inverse  with  elements  v^  »  Then  the  generating  func¬ 
tion  for  these  inverse  elements  is  the  Lagrange  interpolating  poly¬ 
nomial  [30], 


Vs) 


M 

=TT 


(S.8k) 

Jv? 


-1  i-1 


t*  —JL 

■  i  Ti,  • 

i=l 


Sj^Sj  if  k*<J. 


for  which-*1 


V"V  ’  4jk 


Theorem  1  will  assist  us,  in  chapter  III,  in  making  a  direct  com¬ 
parison  of  two  recently  developed  methods  for  exponential  approximation. 


Usually  the  Lagrange  interpolating  polynomial  is  written 

n  ,  v 

TT  8”8k 

Lj(s)  ■  I  j  •  This  difference  is  due  to  the  negative  elements 

k-1  J  k 


{-ai)  in  V 
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III.  DETERMINATION  OF  THE  MATCHED  EXPONENTS  OF  A  DEFINED  ANALYTIC 
FUNCTION 


In  the  previous  chapter,  two  methods  vere  examined  to  find  the 
amplitude  coefficients  in  a  least-square  approximation  of  a  function 
f(t)  by  a  specified  set  of  exponential  functions.  These  linear  least- 
squares  procedures  can  alvays  be  carried  out  given  sufficient  time  and 
precision  to  determine  accurately  the  amplitude  coefficients.  In  con¬ 
trast,  finding  the  complex  frequencies  of  the  set  of  exponential  com¬ 
ponents  to  best  match  the  specified  function  f(t)  is  much  more  difficult. 
Our  attack  on  finding  this  set  of  matched  exponents  begins  with  the 
equations  of  Aigrain  and  Williams  [12], 

III.l  The  Equations  of  Aigrain  and  Williams 

For  a  given  f(t),  t  >_0,  the  necessary  conditions  on  the  2n  para¬ 
meters  to  minimize  the  functional 

i2 


n 


J  *  /  [f(t)  -  l  exp  (skt)]  dt 
o  k=l 

are  expressed  by  the  two  sets  of  n  equations 


(3.1) 


-  /  2  [f(t)  -l  a  exp(s.t)]  [-exp(s  t)]  dt  ■  0, 


3a 


J  o 


k«l 


J 


or 


l  \  I  exp((s  +sk)t)  dt  =  /  f(t)  exp(s  t)  dt 


(3.2a) 


k-1  “  o 


Jsl»2 ,... ,n. 


and 
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-  I  o  /  t  exp((s  +8.  )t)  dt  *•  -  /  f(t)  t  exp(s  t)  dt  (3.2b) 


L  \r  J 

k-1  K  o 


J*l|2 f . . * ,n. 


These  conditions  for  stationarity  of  the  integrated  squared  error  may 


be  written  in  the  frequency  domain 

-l——’  n-B.)  ' 

k-l  V'k  J 


k-1 


2  =  F'("j> 


J-1,2 . n 


(3.3a) 


(3.3b) 


or  even  more  simply  as 

Fa^"8^  “  ^ 

>  F'(.Sj)  J  J-1,2 . n 

where  as  usual* 

00 

F(s)  ■  /  f(t)  expt-st)  dt  (Re{s)>c  ) 
o  ° 

f'(.)  ■  ~  t F( s ) ] 


(3.Ua) 


(3.Ub) 


(3.5a) 


(3.5b) 


F  (s)  ■  j  - 

a  k=l  8“8k 


(3.5c) 


(Equations  (3.1*)  reveal  that  in  approximation  theory,  the  Important 
information  of  the  signal  is  contained  at  the  mirror  images  of  the 

poles  of  F  (s)  which  are  all  points  such  that  Fe{s}>0.  This  suggests 

ft 

that  the  most  useful  information  about  F(s)  and  F  (s)  is  in  the  right 

ft 


SB 

o  »  0  if  /  f^t)  dt 

ft  * 
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half  plane  and  not  at  the  poles  of  F  ( a ) 1  To  further  demonstrate  thlB 

a 

point,  consider  the  following  two  functions, 

g  (t)  ■  exp(-at)  u  n(t)  a  >  0 

g2(t)  ■  exp  (-at)  (u^t)  — u^^( t— t ) J  t  >  0 

Then  the  corresponding  Laplace  transforms  are 


Notice  that  G^s)  does  not  have  any  poles  even  for  arbitrarily  large 
T.  This  means  that  the  pole  of  G^(a)  in  the  left  half  plane  is  due 
solely  to  the  tail  end  of  the  exponential  which  is  a  negligible  part 
of  the  function  g^(t)  for  large  at.)  • 

The  2n  equations  (3.U)  were  formulated  by  Aigrain  and  Williams  in 
19  W.  Despite  their  simple  appearance,  closed-form  solution  of  these 
nonlinear  equations  is  impossible  except  in  trivial  caseB,  Two  ways 
that  have  been  used  to  solve  these  equations  are  gradient  methods  and 
"linear  iterative  schemes".  These  methods  will  now  be  discussed. 

III. 1.1  Gradient  Methods 

One  straightforward  way  of  finding  the  matched  exponents  is  based 
on  the  method  of  steepest  descent.  That  is,  one  finds  a  suitable 
scalar  function  of  the  2n  parameters  {cl^s^}  which  has  a  relative  min¬ 
imum  for  values  of  these  parameters  that  satisfy  the  Aigrain  and  Williams 
equations.  Clearly,  a  suitable  function  is  the  integrated  Bquared-error 
J  defined  in  (3.1).  The  gradient  of  J  is  computed  at  some  initial  point 
in  parameter  space  and  then  the  parameter  point  is  perturbed  in  the 
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direction  of  the  negative  gradient.  The  process  is  repeated  until  the 
gradient  is  approximately  zero. 

Let  f(t)  be  normalized  so  that 

•• 

||f(t)||2  «  <f(t),f(t)>  H  I  t(t)t*(t)  dt  •  1  (3.6) 

o 

Then 

J-  II  f-f  J|2- 1  -2  J  l  \  F»(“8k)*  (3‘7) 

*  *  *  k*l  *  B  * 

The  dependence  of  J  on  the  {a^}  may  be  suppressed  by  using  (3.Ua)  and 

its  equivalent  form  (2.21).  Under  the  constraint  of  equation1,  (3.Ua), 

t  Equation  (3.8)  has  an  interesting  geometric  interpretation.  By 
definition,  the  are  the  coordinates  of  f  on  the  oblique  basis 
B|  whose  elements  are  {expis^t)}.  The  reciprocal  basis  is  defined  as 
I)|  where  |D  ■  |B  (BlB)”1.  A  "vector"  <f  J  nay  be  also  written  as  a 
linear  combination  of  the  dual  basis  elements, 

<*J  "  l  gk  "  <£& 

k=l 

The  square  of  the  length  of  |  is 

IlfJI  2  -  *  <A  B  |  D  G>  -  <A,G> 

n 

"  °k  ®k* 
k*l 

However,  it  is  a  well-known  fact  that  J"||e||  2*  ||  f ||  2-||fJ|  2  in  a 
least-square  approximation.  Hence,  the  {F(-s^) }  are  the  coordinates 
of  f(t)  on  the  reciprocal  basis  of  B|  i.e. 

<fj  ■  X f(-*)  <^[|* 

k»l 
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n 

J  -  1  -  l  a  F(-s  )  (3.8) 

k-1  K  K 

and  from  (2.21) 

J  •  1  -  l  l  F(-s  )  F(-«.)  (3.9) 

k-1  J-l  J 

where  h“j  are  the  elements  of  the  inverse  of  the  Hilbert  matrix  defined 
in  (2.13).  In  matrix  form 

J  »  1  -  F  H-1  F  ■  1  -  F  C  C  F 

-  i  -  (££)  C  F  (3.10) 

where  the  cik  are  defined  in  (2.14).  Hence,  equation  (3.10)  gives  an 
explicit  expression  for  the  integrated  squared  error  in  terms  of  the 
{8k>  only. 

For  the  scalar  function  J,  however ,  gradient  methods  have  two  serious 
shortcomings.  First,  this  error  is  insensitive  to  changes  in  the  {s^} 
and  as  a  result  convergence  to  the  minimum  is  slow.  Second,  because 
of  the  correlation  between  exponentials  the  error  can  be  reduced  to  near 
its  minimum  value  for  a  wide  range  of  {s^},  In  Beveral  cases  tried, 
descent  methods  converged  to  values  other  than  the  minimum.  (Box  [ 31 3 
has  shown  with  several  examples  why  gradient  methods  don't  always  converge 
to  the  minimum.)  As  n  increases,  convergence  to  the  matched  exponents 
by  gradient  methods  becomes  difficult  to  attain  since  the  measure  of 
dependence  between  the  set  of  exponentials  increases  so  rapidly  with  n. 

A  better  method  of  attack,  achieved  by  working  directly  with  the 
Aigrain-Williams  equations,  will  now  be  considered. 
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and  the  {a^}  are  eliminated.  However,  equation  (3.12)  is  hopelessly 
nonlinear  in  (s  }  and  in  its  present  form  has  been  found  to  be  worthless 
for  computing  these  exponents.  The  next  two  theorems  will  help  put  (3*12) 
in  more  useful  form. 

III. 1.2.1  Explicit  Inverse  of  the  Hilbert  Matrix 

In  chapter  II,  a  derivation  for  an  explicit  expression  for  the  in¬ 
verse  of  the  generalized  Hilbert  matrix  was  derived  using  orthonormal 
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exponentials  and  the  restriction  that  the  real  part  of  every  be 
positive.  This  constraint  was  introduced  simply  to  insure  that  the 
exponential  components  could  be  normalized.  Except  for  that,  it  seems 
to  be  unnecessary  and  may  be  removed  by  "analytic  continuation" »  to 
yield  the  more  general  theorem: 

Theorem  2  -  Let  h^^  ■  -l/(s^+Bj)  be  an  element  of  the  n  X  n  generalized 
Hilbert  matrix  H  associated  with  the  set  of  n  scalars  {s^}  with 
a^  Sj  for  all  i  and  J  (i^j)  and  s^O  i«l,2,...,n.  Define  £ 
aa  the  n  X  n  matrix  with  elements 


Us.s . 

■  -A  J-  t"  t11 

V8j  1  * 


where 


5  i  =TT 


k-l 

k?fa 


<VV 

'V.1 


(3.13) 


Then  D  »  H 


-1 


Proof  *2 l  induction 

Let  £  be  the  product  matrix  *  8  D.  We  wish  to  prove  that  ^  is 
the  unit  matrix  with  elements  6^.  Define 


kiJ  2  4tJ  ■  hlk 


?  vx  ^ 

wwv7 ' 


(3.1M 


n-1  .  "j1  ‘V.11T1  ^ 

‘  Ji  Wv7 ' 


From  (3.13)  it  follows  that 
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and  thus 


•t"  .  LVV  t“-1 

k  ^n-k1  k 


.0-1.  f  Vj  ‘‘."i,1  W 

1J  Ji  Tv*?Tvv' 


(3.15) 


n-1  ?  T  T1  f  ^8n“8k^8n“8lOl 

r‘1J  'JiL'v^'vv  l1  ‘  ‘W'VVJJ 


a,n  ?  V.1 

TV“7  Ji  Wvv 


A*1  A0-1 

A . .  -  A.  . 


28. 


ij  u  iv^r 


(3.16) 


For  (3.16)  to  hold  when  J«n  and  i^n,  A^”1  muat  be  defined  to  be  zero. 

For  (3.16)  to  hold  when  i  and  J  are  both  equal  to  n,  one  muat  subtract 
A^1  from  the  right  hand  aide  of  (3.16)  whenever  i*n.  With  this  modi¬ 
fication!  ao  that  it  may  be  applied  generally,  equation  (3.16)  becomes 

*ij  "  Vj1  "  Uin  "  *in  AL1)  i»J“1*2 . n*  (3.l6a) 

n  J  n 


We  now  assert  aa  the  inductive  hypothesis 


“H1- v  (3-i7) 

which  is  readily  shown  to  be  true  for  n«2  and  3.  To  establish  validity 


for  larger  n,  first  substitute  (3.17)  in  (3.16a)  to  obtain 

2s  T11 

Kj-v '  t'A-V  ■  si“1  ■  °-  . -  (3-i8) 

^Sc  *y  n 


Then  for  the  n  X  n  matrix  A 
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AiJ  "  6iJ  +  Un+ijjTn  ^Ain  “  6in* 


(3.19) 


From  the  symmetry  in  equation  (3.14) •  if  the  subscripts  J  and  n  are 
interchanged  equation  (3.19)  must  hold.  Hence, 


4in  '  Sln  *  U1J  '  V 

n  J  J 

Substituting  (3.20)  in  (3.19)  one  gets 
4s  s.  | 

1 - Sj“2  Pij  “  6iJ]  *  ° 

L  'v'j'  J  3  3 


(3.20) 


Thus,  it  is  necessary  that 


■  su 


i »J"1»2 .  ,n 


and  Theorem  2  is  proved  by  induction.  Q*  E«  D« 
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III. 1.2.2  Simplification  of  B  ■  -  G  H  . 

Theorem  3  -  Let  G  and  H*1  be  n  X  n  symmetric  matrices  with  elements 

-48,8, 


e 


13  (Vj> 


and  h~*  -  I  T?  T® 

)2  ij  (s^J  i  J 


y 

-1 


Then  the  elements  of  the  matrix  B  *  —  ^  equation  (3#  15)  are 


b?j--  j  *<>> 


-i 


k=l 

s  T*1 

ll  ^ 


ik  kj 


8i(V9J)T? 


1  *1  , 


Bii"  “  2T  +  l  2 - 2  * 

11  "i  k«l  8,-8, 


(3.21a) 


(3.21b) 


k*i 


’i  k 


Proof 


Formula  (3.21)  holds  for  n*2  since  ■  OH  1  « 


(2a/ 


1 


<v.2!2 


^(•2+81)  |  (2b2)  J 


2s, 


(•jS\ 

■  Va  fi^V 

V  *!-*2  ' 

Vs  /VM*'  /V‘4 

-  Tv^r  (v^j :  2*2(v^) 


1  f  2®2 

82 

281  (•12-«22) 

•Tv*?1 

*1 

i  t  2,i 

•2(«2-Si) 

2*2  fr»Vj 

By  definition 

n 

B?J"  * 

0  k-1 

U  8  T°  T° 

Vj  Ak 

(,i*,k)2(Vj) 

(3.22) 


We  new  establish  an  inductive  principle  to  derive  along  lines  very 
similar  to  that  used  for  Aa  in  equations  (3.1U)  to  (3.16)  it  is 

readily  shewn  that 


un  r11-1 

Bu '  Bu 


B 


New  make  the  inductive  hypothesis 

i-l 


B 


n-1 


8 


in* 


)(8*Ol? 

j-£.  j _ n  y  x 


(3.23) 


lj  ’i(V*jHVjHVi)Ii 

Using  (3.23)  and  (3.2U) 


iyj.  (3.2U) 


Bn 

'  9  T°  1 

J  J 

»  B0"1  ♦ 

'?8  / 

J  J 

BD 

r  st11  l 
Vj 

“ij  * 

BU  * 

L  n  J  nj 

fl.  ■“ 

in 

•i(,i"*j)Ti- 

-s  T° 

JLl± 


1-  >  .  ni  . . 


r  wv^j 


2s,T“  bJ 

1*1  is 

(•  ♦,,)T“ 

n  J  n 


or 


.142. 


r 

V5 

m 

1 

% 

2a  T° 

3.1 

L  •,«•*— j)*?J 

(Vj< 

or 

bh  -  -3  *  A  - 

(b  ♦a.k 
n  J  n 

k- -f*  .1  ■ 

L  *i  *l“*n  ^1-* 


(3.25) 


As  in  Theorem  2,  from  the  symmetry  in  (3.22)  it  is 

B, 


necessary  that 


s  T 


(3.26) 


and  (3.21a)  ia  proved  by  induction.  When  i^J, 

Bn-1  .  n«i  “r1  1  (a  -a.)  n-1  ■  (a  .a  ) 

Bii  ‘  S*i  2  — ■  o  -  -  1*.,  T-S-i.  f  l  *kt  **n  V 

*  a  n  k 


Then 


_  (■-■,)  n  a  T°(a  .a  ) 

-  .  Ua^T?  j-2-X.  V  WW 

V§i  k-1  (•i+"k)3(»n+ek) 

B“  .  B-1  .  .  It*  T?  J  PVVVV  ,] 

11  1 1 A  <w3  •  ll 


(3.27) 


2a. T° 
i  x 


B?  . 


From  (3.26)  and  (3.28) 


B“  -  Bn_1 
ii  ii 


(a  ♦a.)Tn  ln 
n  i  n 


2a 


(3.28) 


,  2  2 

a.  -a 

i  n 


(3.29) 


and  (3.21b)  is  proved  by  induction. 


q.e.d. 
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III. 1.2. 3  Equations  of  the  Direct  Method  in  Integral  Fora 

An  explicit  expression  has  now  been  found  for  the  matrix  B  in 
equations  (3.12)  and  the  equations  of  Aigrain  and  Williams  have  been 
reduced  to  a  set  of  n  equations  and  n  unknowns  in  the  (s^)  alone.  Al¬ 
though  equations  (3.12)  appear  In  the  simplest  form  possible,  they  are 
highly  nonlinear  in  the  {s^}  and  serious  computational  difficulties 
arise  in  their  solution.  With  a  little  manipulation,  these  equations 
may  be  expressed  in  an  integral  form  which,  as  we  shall  see  in  the  next 
section,  has  an  illuminating  geometrical  interpretation.  Even  more  im¬ 
portant,  this  form  is  well  suited  for  solution  using  a  linear  iterative 
method. 

By  multiplying  each  equation  in  (3.12)  by  s^^  one  gets 


Z 

k»l 


s.T® 
1  X 


ik 


F(-V 


0  i«l,2,...,n.  (3.30) 


It  can  immediately  be  shown  that  (3.30)  car  be  written  in  the  more  com¬ 
pact  form 


H( 3 )  F(-s)da 

4-  8-8i  2"J 


0 


i«l,2,...,n 


(3.31) 


where 


n 


H(a)  -  TF 

t=l 


ls+8v) 

(8'sky’ 


(3.32) 


This  is  easily  verified  by  using  residue  calculus  and  noting  that 
„n 


3Tj 


-  !Bii  * 


T 


,n 

r 


(3.33) 


The  relation  of  equation  (3. 31)  to  the  Kautz  procedure  for  construction 
of  orthogonal  functions  and  associated  development  by  others,  is  discussed 


next. 
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iii. i.3  Method  of  McDonough  and  Huggins 

In  a  recent  paper  [7]  McDonough  and  Huggins  suppressed  the  ampli¬ 
tude  coefficients  (o^)  in  the  Aigrain-Williams  equations  by  regarding 
f(t)  as  a  signal  in  a  vector  space.  Their  argument  proceeded  as  follows: 

Let  the  error  of  the  approximation  be  e(t)®f(t)-f  (t).  Then  it  is  read- 
ily  seen  that  equations  (3.2a)  and  (3.2b)  may  be  written  respectively  as 

OB 

J  e(t)  expfs^t)  dt  ■  0 
o 

and 

QO 

J  e(t)  t  exp(skt)  dt  ■  0 
o 

In  the  language  of  vector  spaces,  equation  (3.31*)  means  that  the  error 
e(t)  is  orthogonal  to  the  space  S2n  which  is  spanned  by  the  2n  "vectors” 
(exp(s^t),  t  exp(s^t) }.  Also,  by  definition,  the  approximating  func¬ 
tion  f&(t)  must  lie  in  the  subspace  Sr  spanned  by  the  n  vectors 
(exp(sjt)}.  Let  s2n“Sr.  ^encrte  *he  oubspace  of  S2r  that  is  complementary 
to  Sn  and  let  ^(t)}  i"l,2,...,n  be  basis  for  this  subspace.  A  basis 
for  S2n"sn  can  be  formed  by  applying  the  Gram-Schmidt  procedure  to  the 
functions  (ax^'s^) ,. . .  ,exp(snt) ,  t  exp^t) ,. . .  ,t  exp^t)  )  taken  in 
that  order,  to  construct  the  orthonormal  basis  functions  i“l,2,. .. ,2n. 

Clearly.  {$  }  i»l,2,...,n  is  orthogonal  to  both  f  (t)  and  e(t)  and 

n+i  ...  a 

hence,  must  be  orthogonal  to  f(t)  *  e(t)  +  f&(t).  The  $n+i(s),  which 
are  the  Laplace  transforms  of  <f>n_^(t),  may  constructed  by  simple  in¬ 
spection  from  Kautz's  method.  That  is, 


$  ^.(s)  =  H(s)  $.(s) 
n+i  l 

where  again 

H(s)  * 

k=l  '“-“k 


JJ  (8+8k) 


(3.35) 


(3.36) 
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and 


•cV’i  TTlflid 
(■-O 


♦1(s) 


•-B, 


(3.37) 


k-1 


Now,  n  independent  equations  of  constraint  on  the  (s.  }  must  be 


J  Vi(t)  dt 


i*l ,2 | • • . ,n, 


(3.38) 


which  when  written  in  the  frequency  domain,  using  the  Parseval  relation, 
are 

Jo. 

i-1,2 ,. . . ,n.  (3.39) 


J* 

I 

-J 


L  F(-‘)  *n*l(,)  Hls)  5Sj  ”  ° 


Since  these  equations  involve  the  {s^}  only,  the  {a^}  have  been  suppressed 
as  in  (3.31).  In  fact,  it  will  be  shown  in  the  next  section  that  the 
left-hand  sides  of  equations  (3.39)  are  merely  linear  combinations  of 
the  left-hand  sides  of  equations  (3.31).  These  equations,  (3.31)  or 
(3. 39) »  are  still  nonlinear  in  terms  of  the  unknowns  (s^}  and  apparently 
one  of  the  b-*8t  ways  to  solve  for  them  is  by  a  numerical  linear  itera¬ 
tive  scheme  first  suggested  by  Sears  [32]  and  described  as  follows: 

Let  the  all-pass+  operator  H(s)  in  (3.39)  (or  (3.31)  )  be  replaced 


by  the  more  general  operator 

n+1  b,  s*'1  (-l)nD  (-s) 

- TO — 


H*(8>  ’  j i  W 


,  b  -  i 
’  n+1 


(3.Uo) 


where 


D(s)  -  (s-sk). 


(3.1*1) 


k-1 


H(s)  is  sometimes  called  the  all-pass  operators  for  if  s  is  replaced 
by  Ju.  the  magnitude  of  H(s)  is  one,  independent  of  the  value  of  the 


frequency  u 
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If  H  (s)  is  used  in  (3.39)  instead  of  K(s),  one  obtains  n  simultaneous 
a  * 

equations  which,  being  linear  in  the  (b^,b2,...  .b^},  may  be  written 
in  matrix  form  as 


M  B  -  Z 


where  M  is  an  n  X  n  matrix  with  i,x  element 
mik  -  /  [F(s)  eW/D(.)l  ^(s)  ^ 


ds 

2nJ  * 


(3.1*2) 


(3.1*3) 


B  and  Z  are  columns  with  elements 


Zi  “  "  “i,n+l 


(3.UU) 


and  are  the  unknown  coefficients.  The  iterative  algorithm  is  as 
follows : 

(a)  Given  the  poles  at  the  J*'*1  iteration,  i.e., 

{s1,s2,...,sn}j 

evaluation  the  matrix  (M),  and  the  vector  (Z).. 

=  J  ”  J 

(b)  Solve  equation  (3.1*2),  to  obtain  the  coefficients  of  the 
vector  (B) 

(c)  From  (b)j+1  find  the  new  pole  locations  (s^,, . , ,sn>j+1  using 

(d)  Repeat  from  (a)  with  J=j+1.  Continue  the  process  until 

the  change  mfX| (s. ),-(s, ) .  ,1  is  less  than  some  small 
l  l  j  l  j+i 

pre-assigned  value. 

The  convergence  properties  of  this  method  shall  be  discussed  in  chapter 


III.  1.1*  Equivalence  of  the  Direct  Method  and  McDonough's  Method 


Hf.  (s)  -  H(s)  -±- 
1  s-s. 


i  =  l  |2  |  •  •  •  |Zs  • 


(3.1*5) 


Then  fran  Kautz's  method  it  is  easily  seen  that  the  {^(s))  also 
form  a  basis  for  the  difference  space  Sgn-Sn  and  hence ,  equations 
(3.31)  have  the  same  geometric  interpretation  used  by  McDonough. 

Thus, 

00 

/  F(-s)  ^(a)  =  o  i-l,2,...,n  (3.46) 

o  " 

are  Just  linear  combinations  of  equations  (3.39).  Notice  that  the 
all-pass  function  H(s)  is  still  preserved  and  the  same  linear  itera>- 
tive  scheme  can  be  used. 

These  new  equations  in  terms  of  ¥^(s)  have  two  advantages  over 
(3.39).  First,  ^(s)  has  only  one  double  pole  whereas  *n+i(s)  has  i 
double  poles.  This  means  the  old  set  will  have  i-1  extra  derivative 
terms  when  the  residues  are  evaluated  and,  moreover,  each  term  will 
have  (i-1)  extra  factors  of  the  form  Clearly,  there 

is  a  saving  in  computational  time  by  using  the  new  equations.  Second, 
equation  (3.46)  may  be  written  in  the  matrix  form 

F'»BF.  (3. 47) 

This  enables  one  to  use  matrix  algebra  to  find  the  (s^)  using  (3.21). 
Solution  of  the  equation  in  this  form  has  not  been  attempted  here,  but 
is  a  topic  for  further  investigation. 

While  the  *n+i(a)  orthonormal,  the  Y^(s )  are  not.  As  a  result J 

+  One  is  contrasting  an  orthogonal  versus  an  oblique  basis,  both  of  which 
span  the  same  space.  Naturally,  there  will  always  be  more  correlation 
between  the  oblique  elements.  However,  here,  as  is  usually  the  cose, 
the  oblique  elements  are  easier  to  express  mathematically  and  any  gains 
in  accuracy  made  by  using  orthonormal  elements  may  be  lost  due  to  the 


extra  complexity  introduced  into  these  expressions 


— U8— 


although  equations  (3 .1*6)  have  the  simpler  form 


n  r  J°°  k-1  ,  ,  1  J“  _n  . 

L  Li  n's)  3lJ  V;/  F(-‘>  dI'DI.:'.") 


k-1  I--J 
n 


or 


Jx 


P  B  ■  R, 


i**l,2,...,n 


(3.1*8) 


they  are  "softer"  than  equations  (3.1*2). 

Evaluating  the  elements  of  the  matrix  P  by  residue  calculuB ,  one 

gets 


( 


'ik 


-’i1"'1  X  T^prT}F(-8i>. 

m^i 


where 


n  (-s)*"1  P(-s  ) 
*  T  m  m» 

L  (a  -8.)  r 

m=l  m  l  m 

m^i 


it 

r.  ■  TT  (s  -a  ) 

1  m=l  1  B 
m^i 


(3.1*9) 


(3.50) 


Since  f(t)  is  real,  the  {s^}  must  occur  in  complex  conjugate  pairs. 
Upon  examining  (3.1*9)  it  is  seen  that  if  s^  is  replaced  by  s*,  p^  becomes 
p^.  ThiB  also  implies  that  if  Si  is  real,  so  is  p^.  Hence,  it  ir 
possible  to  avoid  complex  arithmetic  altogether  in  finding  the  real 


{bi>  from  (3.1*8)  by  using  the  equivalent  system  of  n  equations 


l  pu  ■  ri 


J-l 


i=l,2,.. . ,NREAL 


(3.51*) 


and  Williams.  The  key  feature  of  the  MSS  method  is  the  introduction 


of  an  approximate  error  E  (s),  viz. 

a 


El. (a)  N  (a) 

V>-3^rr«->-5^r 


(3.5U) 


where  now  the  subscript  j  refers  to  the  iteration  number.  To  solve 

the  equations  in  a  feasible  way,  the  previously  computed  (b,  )  co- 

K  J  *  Ji 

efficients  of  D  ^  are  regarded  as  fixed  during  the  J  iteration. 

This  linearizes  the  error  in  terms  of  the  unknown  coefficients 

{a  ,b  }  of  the  numerator  polynomials  N  and  D,.  One  now  replaces 
K  K  J  J  J 

e(t)  in  (3.53)  by  e  (t),  the  inverse  transform  of  E  (3),  and  uBea  an 

a  a 

iterative  scheme  very  similar  to  the  one  employed  by  McDonough  described 
earlier.  With  repeated  iterations  D^Cs)  approaches  1)  (s)  and  thus 
L  (s)  approaches  the  true  error  K(s)«F(s)-F  (s). 

8-  A 

However,  inserting  e  (t)  in  (3.53)  has  three  distinct  diB advantages. 

A 

First,  instead  of  utilizing  the  convenient  point  form  of  the  Aigrain- 
Williams  equations,  one  must  evaluate  a  set  of  2n  partial  derivatives 
and  then  integrate.  The  resulting  equations  are  much  more  complicated 
than  (3.M.  Second,  and  more  important,  the  Mode-1  Iteration  used  on 
these  2n  new  equations  does  not  converge  to  the  optimum  solution  since 
the  approximate  error  is  minimized  rather  than  the  true  error.  Thus, 
following  the  Mode-1  Iteration,  a  Mode-2  Iteration  is  sIbo  needed  to 
further  refine  the  results  and  find  the  optimum  solution.  This  diffi¬ 
culty  does  not  arise  in  the  Mode-2  Iteration  bece.ose  the  expressions 
are  correct  in  the  limit  as  the  (bv)  approach  their  optimum  values. 
Third,  several  examples  will  demonstrate  that  Mode-2  converges  more 
slowly  than  McDonough's  method  and  a  new  method  which  we  will  present 


i"l  ,2, • • • ,n. 


Hence,  equations  (3.58)  provide  another  linear  iterative  scheme  involv¬ 
ing  2n  real  parameters  (a^,^)  instead  of  Just  the  n  {bk>.  Otherwise, 
the  iterations  are  carried  out  as  in  the  McDonough  method. 
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Mode-2  Iteration  of  the  MSS  Method 
The  error  E(s)  may  be  written 

E(s)  ■  F(s)  -  F&(s)  *  F(a)  -  . 


(3.59) 


Thus 


islsi.  4^ 

3a.  D(s) 


3E(s)  _  rg^Kjs)  -s1"1  _  (  v 
D2(a)  a 


3b. 

l 


(3.60) 


Using  the  Parseval  relation  on  (3.53)  gives  the  conditions  for  the  sta- 

tionarity  of  the  integrated  square  error  in  the  frequency  domain  as 
J“  -  «./  _  \*  •»  _i~l 


-J 

J 


i  [F(-s)  -fe}][§rrr]  rff' 0 


I  CF(-8)  -  m  t 


i-i 


-21. 

0 

(3.61a) 

2irJ 

§$■] 

£|& 

B 

O 

(3.61b) 

i*l (2  § • • • 

If  the  iterative  process  for  Mode-1  converges,  Dj_^(s)  approaches  Dj(a) 
and  comparison  with  (3.6la)  shows  that  (3.58a)  is  correct  in  the  limit. 
However,  (3.58b)  is  not  correct  in  the  limit  which  ia  observed  when 
it  is  compared  with  (3.6 lb ) •  For  this  reason,  Mode-1  Iteration  does 
not  converge  in  general  to  the  optimum  solution.  This  difficulty  may 
be  eliminated  by  using  (3.60)  to  change  (3.58b)  to 


f  r°,,<->F<-°>-V">|  r_£L.  W! 

4-  L  Vi(s)  J  L vT71  vT* 


i-1  N . 


,(•)! 
17 


da 

2wj 


(3.62) 


Equations  (3.58a)  and  (3.62)  are  now  used  in  the  Mode-2  Iteration. 
Convergence  to  the  optimum  solution  is  now  usually  possible,  but  as 
we  shall  see  in  chapter  IV,  Mode-2  converges  so  slcwly  that  in  order 
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to  make  the  MSS  method  practical,  one  must  first  use  the  more  rapidly 
converging  Mode-1  Iteration  to  bring  one  "near  enough"  to  the  optimum 
point  in  parameter  space.  Furthermore,  there  is  no  guarantee  Mode  2 
converges  if  one  is  not  "near  enough"  since  the  equations  that  de- 
termine  it  are  not  correct  unless  one  ia  actually  at  the  minimum. 


III. 1.6  The  Nev  Method 

Thus  far,  ve  have  discussed  two  linear  iterative  schemes  in 
sections  III. 1.3  and  III. 1.5.  Each  has  worked  well  for  the  cases 
reported  and  appears  to  be  useful  in  finding  by  clerical  computation 
the  matched  exponents  for  the  approximation  of  a  known  time  function. 
Ibis  section  develops  yet  another  linear  Iterative  method  which  offers 
the  advantages  of  both  the  methods  described  in  sections  III. 1.3 
and  III. 1.5  and  reveals  the  link  between  them.  Ibis  method  leads 
to  the  same  results  as  those  described  by  McDonough  and  Huggins. 


idamental  Equations 


The  Aigrain-Williams  equations  (3.k)  may  be  written  in  the  form. 


E(-.k)  -  F(-S.t)  -  Fa(-sk)  «  G 
E'(-sk)  -  F'(-sk)  -  F'(-«k)  «  0 


(3.63) 


The  equations  in  this  form  suggest  that  a  better  way  of  using  the  ap¬ 
proximate  error  ^&(s)  defined  by  (3.5*)  is  to  impose  che  constraints 

of  equations  (3.63)  directly  upon  it.  This  lasted! ately  yields  a  set 

th 

of  new  linear  equations  for  the  J  '  Iteration. 


E^(a)  ■  0 

e;(s)  -  o 


for  s  »  (—8,  ) 

k  4-1 


k»l ,2,. . . ,n 


(3.6km) 

(3.6kb) 


obtained  from 
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where  the  are  the  roots  of  the  denominator 

the  previous  iteration.  Equations  (3.6Ua)  may  be  expressed  as 

D  F-N  =  0,  s  *  (-s  )  k*l,2,...,n  (3.65a) 

Similarly,  upon  differentiating  E  (s)  with  respect  to  s,  equations 

ft 

(3.6to>)  may  be  written  as 

E'-  (DJ.1(D/%DJ»  . 

-  0,  »  -  k-1.2 . ” 

which  by  using  (3.65a),  simplifies  to 

r'YF  Dj  ■  "j  ■  5  ■  (-"k >>1  k’1*2 . n  (3-65,>’ 

The  2n  simultaneous  equations  (3.65)  are  linear  in  terms  of  the  unknowns 
(a^,b^}  .  The  iterative  procedure  is  carried  out  in  the  same  way  as 
described  in  the  tvo  previous  methods  using  equations  (3.65)  or  the 
equivalent  equations  (3.67)  for  convenient  computation  to  find  the 
{ak  ,bk}  .  The  initial  poirt  in  parameter  space  may  be  determined 

J 

perhaps  from  Prony’s  method  or  Pade  approximants. 

A  very  important  consequence  of  imposing  these  constraints  upon 
E  (a)  is  that  the  D  .  required  in  the  formulation  of  the  MSS  method 
does  not  appear  in  equations  ( 3.65) ;  the  introducti  m  of  the  approvi* 
mate  error  Ea(s)  was  unnecessary.  In  fact,  an  appropriate  set  of 
linear  equation  may  be  gotten  directly  from  the  Aigrain^Wi Ilians 
constraints  as  follows: 

From  (3. ‘•a)  and  (3.52) 

N(s)  *  D(s)  F(s)  for  s  ■  (-s^)  i»l,2,...,n  (3.66a) 

where  the  (s^ )  are  the  roots  of  the  n*h  degree  polynomial  D(s).  Also, 

F'  •i-D*  N  ♦  »'  D)/D2 
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Frcoi  (3.68),  (3.69)  and  Theorem  U,  £  can  be  written  explicitly  in 


terms  of  the  { s ^ }  as 


;  -  h  =  V  S  (—a  ^  W_a  \ 

ik  ik  i  Hij  i'  "  _i' 


(3-70) 


and  U  may  be  regarded  as  the  negative  of  the  n+1  n  column  of  £  or 


u  ■  c 
i  i  ,n+l 


i"l,2 | . . . ,n. 


(3.71) 


But  (3.U8),  (3.!»9),  and  (3.70)  reveal  that 


ri  fij  ■  cu 

and  thua  McDonough's  method  and  the  one  developed  here  must  be  equivalent. 
Discussion 

In  this  section  we  have  revealed  the  strong  link  between  the  MSS 
method  and  that  of  McDonough.  Both  methods  use  equation  (3.5*0  or  its 
equivalent  to  linearize  the  iterative  process  (although  this  was  not 
so  obvious  in  the  latter).  The  crucial  difference  in  the  methods  is 
that  the  MSS  method  considers  variation  of  the  error  with  respect  to 
the  ( a^,b^}  parameters,  whereas  in  McDonough's  method  (and  the  one 
developed  here),  the  variation  with  respect  to  the  exponents  {s^}  (and 
the  hidden  (a,  })  is  considered.  It  is  not  possible  to  write  a  set  of 

it 

linear  equations  in  the  {a^b^}  for  the  true  error  surface  whereas 
equations  (3.65)  and  (3.66)  show  that  one  may  do  this  when  the  varia¬ 
tion  is  with  respect  to  the  {0^,8^},  thus  avoiding  the  need  for  two 
types  of  iterations. 


In  conclusion,  the  new  method  developed  here  has  several  computa¬ 
tional  advantages  over  the  one  developed  by  MSS.  First,  it  requires 
only  one  iterative  scheme  instead  of  two.  Second,  by  using  the  point 
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Torm  of  the  Aigrain-Williams  equations,  the  matrices  in  equation  (3,67) 
immediately  appear  as  explicit  function  of  the  {•  }-  Tn  the  MSS  method 
the  corresponding  matrix  elements  (see  Table  U.i,  p.  6l)  are  much 
harder  to  evaluate.  Thirl,  as  we  shall  show  in  chapter  IV,  for  all 
examples  thus  far  examined,  the  new  method  converges  more  quickly  to 
the  matched  exponents  than  the  MSS  method. 
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IV.  CONVERGENCE  AND  COMPARISON  OF  THE  LINEAR  ITERATIVE  SCHEMES 
In  the  previous  chapter  three  linear  iterative  schemes  vere 
described.  Two  of  these,  the  method  of  McDonough  and  the  new  method, 
vere  shown  to  be  equivalent  in  their  results,  although  computationally 
different  I  That  is,  for  any  initial  {s^}  and  fixed  n,  the  resulting 
Iterations  of  either  of  these  methods  will  be  Identical  barring 
roundoff  errors.  However,  because  the  new  method  uses  the  rational 
fraction  form  of  F  (s),  a  direct  comparison  of  it,  rather  than  McDon- 
ough's  scheme,  with  the  MSS  method  will  be  made  since  this  will  be 
much  easier  to  do. 

IV. 1  Comparison  of  the  Iterative  Equations 

The  2n  equations  used  in  the  iterative  scheme  of  MSS  may  be 
written  in  matrix  form  as 
VI  A  ♦  01  B  ■  XI 


W1  A  +  HI  B  -  Y1 


U.l) 


for  Mode-1  Iteration  and  as 


VI  A  ♦  G1  B  -  XI 


W2  A  ♦  H2  B  ■  Y2 


(U.2) 


for  Mode-2  Iteration.  From  (3.56a)  it  is  seen  that  the  elements  of 
VI  are 


l-.)*-1  B1’1 


“  " ' ■’jJ-i'vim  ST 


,n-*  n  MO 


vi-Hs-3 


TT 


(k.3.) 


i |K9jl i^|»«  •  §n j 


compared  with  the  corresponding  ouch  simpler  expression  v^«(-s^r 

given  in  equation  (3.67)  for  the  new  method,  and 

-1  .  f°°  sk“1F(s)  ds 

ik  L  D. 


-S-  Vi(-8)  Vi(8)  2*J 


xli  “  “gli ,n+l* 

From  (3. 5$b )  it  is  shown  that  the  remaining  elements  in  equation 


(4.3b) 


(U.3c) 


(4.1) 


ik  L  D,  ,  -s) D,  ,(s)  2itj 


-J~  J-l 


(— s  s^“^F(b)  F(-s)  ds 


hi  -  -  f  lr-sl  ■--? — £j&i- 
“  -J-  Vi(s) 


(4.4a) 


(4.4b) 


yli  "  “  hi  ,n+l* 


(\.hc) 


Equation  (4.4b)  can  be  difficult  to  evaluate  by  residue  calculus. 
For  example,  if  f(t)  is  the  square  pulse,  F(s)  *  (l-e""8)/s,  the  product 
of  F(s)  F(— a )  with  any  rational  function  of  s  will  have  an  essential 
singularity  at  infinity,  in  both  the  right  and  left  hand  planes,  and 
thus,  direct  evaluation  of  (4,4b)  by  residues  requires  special  treat- 
mentj  (Of  course,  one  may  evaluate  these  Fourier  transforms  by  direct 
integration  with  respect  to  w  over  -•  to  “  without  resorting  to  resi¬ 
dues  but  this  is  usually  extremely  difficult.)  Alternatively,  one 
may  invoke  the  Parseval  relation  and  evaluate  the  equivalent  time-domain 
equations  (3.56)  to  obtain  the  matrix  elements  for  Mode-1  Iteration. 

Table  4.1  summarizes  all  of  these  equations  and  clearly  shows  that 
the  elements  in  the  new  method  are  much  easier  to  compute. 


t 


See  Appendix  C  for  details 


Table  U.l  -  Comparison  of  the  Methods 


r 
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IV. 2  Rates  of  Convergence 

To  predict  by  mathematical  analysis  the  rate  and  region  of  con¬ 
vergence  of  any  of  the  linear  iterative  schemes  for  a  general  f(t) 
is  extremely  difficult  except  for  the  simple  case  when  n*l.  Instead, 
ve  provide  severed  numerical  examples  to  give  the  reader  seme  feel 
for  results  obtained  by  the  different  methods.  For  any  f(t)  that 
is  composed  exactly  of  n  exponentials,  any  of  the  iterative  processes, 
Mode-1,  Mode-2,  and  the  new  method,  will  yield  these  exponentials 
immediately  after  one  iteration.  Consequently,  when  f(t)  is  "nearly" 
exponential,  one  would  also  expect  reasonably  rapid  convergence  for 
any  of  the  methods.  This  is  indeed  the  case  as  ve  now  illustrate  by 
specific  examples. 

numerical  Results  -  Consider  the  two  time  functions 

f^t)  *»  (e-t  +  e“2t)  u^t)  (4. 5) 

and 

f2(t)  «*  (e_t  -e”2*)  u^(t)  (4.6) 

each  to  be  approximated  by  a  single  exponential.  Figure  1  shov^  quali¬ 
tatively  why  f^(t)  can  be  approximated  very  accurately  by  one  exponen¬ 
tial  whereas  f2(t)  cannot. 

Table  4.2  gives  the  result  for  f^(t)  of  the  iterations  by  the  var¬ 
ious  methods  all  starting  with  on  initial  value  of  s^*  -  1.2.  The  new 
method  and  Mode-2  Iteration  converge  to  the  same  result  (the  optimum 
approximation )  but  Mode-2  required  >U0  iterations  whereas  the  new  method 
needed  only  4.  In  contrast,  Mode-1  converged  as  rapidly  ao  the  new  method, 
but  not,  to  the  optimum  approximation.  Bjus,  for  this  simple  signal,  the 
new  method  is  superior  to  both  Mode-1  and  Mode-2. 
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Flgure  1  —  Plot*  of  f^ (t)  and  f2(t)  and  the  heat  approximation*  of  them 
by  a  single  exponential.  d 


For  f2(t),  shown  in  Table  U.3,  the  Mode-2  Iteration  took  about 
3000  cycles  I  Also,  observe  that  the  Mode-1  error  is  considerably 
larger  than  in  the  case  for  f^(t).  Ibis  ia  reasonable  since  f  (t) 
more  "nearly"  resemble*  a  single  exponential  than  uoes  f 2 ( t ) .  (Recall 
that  Mode-1  only  gives  exact  values  when  f(t)  is  an  exponential.) 

From  equation  (3.51*),  the  linearized  error  for  approximating  the 
square  pulse,  F(sMl-e“*)/s,  by  a  single  exponential  becomes 


(bJ.  -(b)  t 

e.lt) '  T*Tv7  ' 


[(b  )  -(b  )  ]  -(b  )  t  -(b  ) 

- V  y  1  >-■  e  1  -5"1  (u  .(t)-e  1  ^  u  .(t-l)} 

v  rj-i 

(M) 


where 


f^(t)  ■  u_1(t;-u<_1(t-i). 

Table  U.U  reveals  that  Mode-2  Iteration  and  the  new  method  converge  to 
the  optimvm  exponent  at  the  same  rate.  (By  coincidence,  the  iterations 

are  almost  identical  for  this  oase.  They  diffe.  in  8th  or  9th  decimal 
place,) 


Table  4.2  -  Approximation  of  the  function 
f1lt)=lexp(-t)+exp(-2t)]u<iil(t)  by  One  Exponential. 


J 

‘Vj 

(Vj 

(*i>j 

<Vj 

‘1’j 

<Vj 

1.20000 

1.20000 

1 .20000 

1 

1.84197 

1.20139 

1.93499 

1.32265 

1.93369 

1.32095 

2 

1.84309 

1.20277 

1.93779 

1.32639 

1.93908 

1.32815 

3 

1.8442 

1.20416 

1.93787 

1.3265 

1.93938 

1.32856 

4 

1.84531 

1.20555 

1.93788 

1.3265 

1.9394 

1.32859 

5 

« 

% 

• 

• 

• 

1.93788 

1.3265 

1.9394 

1.32859 

100 

• 

1.92123 

« 

1.30392 

101 

1.9216 

1,30442 

102 

1.92197 

1.30492 

103 

1.92233 

1, 30541 

104 

• 

1.92269 

• 

1.30588 

• 

• 

• 

200 

• 

• 

1.93739 

• 

1.32584 

201 

1.93744 

1.3259 

202 

1.93748 

1.32596 

203 

1.93753 

1.32602 

204 

a 

1.93757 

a 

• 

1.32608 

t 

t 

• 

1.93919 

• 

1.32831 

sE  1  wofiil 

1.9392 

1.32831 

1.9392 

1.32332 

1.93921 

1.32833 

eE 

1.93921 

» 

1.32833 

a 

a 

400 

9 

1.93938 

• 

1.32856 

401 

1.93938 

1.32856 

402 

1.93938 

1.32856 

403 

• 

0 

1.93938 

9 

• 

1.32856 

a 

a 

0 

500 

• 

1.9394 

• 

1.32859 

Mode-2  Iteration 

Mode-1 

Iteration 

Nev 

Method 
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Xable  1,3  -  Approximation  of  the  function 
f1(t)*[exp(-t)-exp(-2t)]u_1(t)  by  One  Exponential. 
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In  this  example  the  term  containing  the  facter  [(b^) j— ("b^ ) has 

little  affect  on  the  equations  that  determine  the  iterations.  However, 

it  seems  quite  possible  that  the  extraneous  terms,  which  always  arise 

when  using  Mode-1  or  Mode-2,  that  contain  the  factor  [(b^)j-(b^) 

could  sometimes  affect  the  equations  enough  to  prevent  convergence  of 

Mode-2  if  (b.)  is  not  "near  enough'*  to  its  optimum  value. 

^  J 


Table  4.4  -  Approximat  ion  of  the  Square  Pulse  by  One  Exponential 


J 

(alJJ 

<Vj 

(v. 

(Vj 

(\>j 

(Vj 

1.00000 

1.00000 

1.00000 

1 

1.5121? 

1.39221 

1.38344 

1.1885? 

1.51217 

1.39221 

2 

1.39272 

1.188 

1.3638 

1.14261 

1.39272 

1.138 

3 

I.4511 

1.29183 

1,36851 

1.15348 

I.45II 

1.29183 

4 

1.42045 

1.23836 

1.36739 

1.15089 

1.42045 

1.23836 

5 

I.43597 

1.26572 

1.36766 

1.15151 

1.43597 

1.26572 

6 

1.42796 

1.25167 

1.36759 

1.15136 

1.42796 

1.25167 

7 

1.43206 

1.25887 

1.36761 

1.15-4 

1.43206 

1.25867 

8 

1.42995 

1.25518 

1  3676 

1.15139 

1.42995 

1.25518 

9 

1.43103 

1.25707 

1.36761 

1.15139 

1.43303 

1.85707 

10 

1.43048 

1.2561 

1.36761 

1.15139 

1.43048 

1.2561 

11 

1.43076 

1.2566 

1.43076 

1.2566 

12 

1.43061 

1.2563** 

1.43061 

1.25634 

13 

1.43069 

1.25648 

1.4-W69 

1.25648 

14 

1.43063 

1.25641 

1.43065 

1.25641 

15 

1.43067 

1.25644 

1.43067 

1.25644 

16 

1.43066 

1.25643 

1.43066 

1.25643 

17 

1.43067 

1.25643 

1.43067 

1.25643 

18 

,  1.43066 

1.25643 

1.43066 

1.25643 

19 

1.43066 

1.25643 

1.43066 

1.25643 

Mode-2  Iteration  Mode-1  Iteration  New  Method 


Table  4.5  shows  the  results  fitting  a  square  pulse  using  3  exponen¬ 
tials  with  the  initial  values  of  the  parameters  chosen  as  (i.)  * 
{-l,-2,-3K  After  90  iterations  the  new  method  had  converged  to 


s  ■  -2.21*6602,  $  ■  -1.1*1*361*3  ♦^4,150741  which  is  in  agreement  to 

i  2 ,  j  ~ 

6  significant  figures  with  McDonough's  result,  ll]  p.  159,  found  by 
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a  search  method. 


Table  1*«5  -  Three  Exponential  Approximation  of  the  Square  Pulse. 
Iteration  Exponents  After  That  Iteration 


-3.0 

-2.0,  -1.0 

1 

-3.01*7813 

-2.037829  ♦  J3.U33861* 

2 

-2.065591 

-1.153318  ♦  J3. 881929 

3 

-2.5U0835 

-1.7860U2  “  jl*. 0946U8 

1* 

-2.050865 

-1. 195975  “  Jl*. 093363 

5 

-2.1*37351* 

-1.686921  ♦  Jl*.  160736 

6 

-2.091982 

-1.21*6309  “  JU.125865 

7 

-2.385592 

-1.62531*9  ♦  Jl*.  165192 

8 

-2.129335 

-1.296223  ♦  Jl*.  131*689 

9 

-2.350106 

-1.580867  ♦  Jl*.  162651 

10 

-2.15831*6 

-1.33U265  ♦  Jl*.  139171 

11 

-2.32U005 

-1.5U77U1  ♦  Jl*.  159901* 

12 

-2.180300 

-1.362868  ♦  Jl*.ll*2197 

13 

-2.301*550 

-1.522938  “  Jl*.15769l* 

Ik 

• 

-2.196827 

• 

-1.381*329  ♦  Jl*. 11*1*391 

• 

• 

k5 

• 

• 

-2.21*7180 

• 

-1.1*1*9385  ♦  Jl*.  150813 

1*6 

-2.21*6105 

-1.1*1*8001  7  jl*.  150679 

1*7 

-2.21*7036 

-1.1*1*9200  7  Jl*.  150795 

1*8 

-2.21*6230 

-1.1*1*8162  7  Jl*.  15069I* 

1*9 

-2.2U6927 

-I.UU9060  7  Jl*.  150781 

50 

• 

-2. 2U6 323 

e 

-1.1*1*8283  7  Jl*.  150706 

a 

• 

82 

e 

• 

-2.2U6600 

• 

-1.1*1*8639  ♦  J»*. 1507l*l 

83 

-2.2U6607 

-1.1*1*861*7  7  Jl*.  15071*2 

8U 

-2.21*6602 

-1.1*1*861*0  ♦  Jl*.  1507U1 

IV. 3  Accelerated  Convergence  -  Shanks*  Method 

Although  the  new  method  converged  to  the  optimum  solution  faster 
then  the  MSS  method  (even  assuming  a  svitch  from  Mode-1  to  Mode-2  is 
made  at  the  best  iime,  a  decision  which  is  apparently  ad  hoc)  in  every 
case  tested,  it  too  converged  rather  slowly  in  some  cases  -  90  iterations 
with  n«3  for  the  square  pulse  and  2U  iterations  with  n«l  for  fg(t).  For 
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larger  n  it  aeons  that  convergence  would  be  e  en  s lower.  In  an  attempt 
to  speed  up  convergence  one  may  incorporate  a  little  used  method  due  to 
Shanks  [ 35 ] . 

Consider  any  of  the  numerical  sequences  in  Tables  U,2  through  U.5 
which  are  either  monotonic  or  oscillatory.  Draw  a  smooth  curve  through 
these  discrete  points.  Typical  graphs  are  depicted  in  figure  2. 

Figure  2  —  Graphs  demonstrating  transient  characteristics  in  the  itera¬ 
tive  sequences. 


New  Method  —  f^t),  f  (t)  New  Method  —  f^t) 

Mode-1  —  f2(t),  f3(t)  Mode-1  —  f^t) 

Mode-2  —  f  (t)  Mode-2  —  f^t),  fg(t) 

These  graphs  look  like  first- or  second-order  transients.  By  "nth  order 

transient"  we  mean  any  function  which  has  the  form 

n 

p(t)  *  B  ♦  l  c,  exp(«a  t)  Re  (a  }  >  0. 

1=1 

Shanks'  method  predicts  the  limit  B  of  such  sequences  by  "filtering  out" 
or  annihilating  the  exponential  components. 

Tables  U.6  shews  the  result  of  applying  Shanks'  method  to  2  sequences 
obtained  by  using  the  new  method.  In  both  cases  the  thirteenth  iteration 
is  correct  to  only  two  decimal  digits  or  so,  but  is  already  correct  to 
six  digits.  In  any  event,  extreme  caution  must  be  exercised  when  applying 
Shanks'  method  to  these  sequences  since  there  is  no  theory  to  Justify  its 
use  here.  However,  it  has  been  demonstrated  how  helpful  the  method  can 
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a  onetime a  be  in  reducing  the  number  of  iterations  needed  and  is  a  topic 
worthy  of  further  investigation. 

Table  U,6a  -  Shanks'  Method  Applied  to  the  First  13  Iterations  of  the 
Matched  Exponent  of  a  One  Exponential  Approximation  of  f2(t)  by  the 
New  Method. 


J 

<v< 

V 

el 

C2 

*3 

% 

*5 

\ 

5.ooooonou  oo 

6.9?  30 /Mil- 01 

9.?l4993t(i-ol 

4. 56290240- 01 

4.574  70990-01 

4.574  5^290-01 

7 

1.  749?  J001.  00 

ft.fl**  JO^/HU-ni 

4.32322 /4D-01 

4.572594)0-01 

4.57429650-01 

4.57430170-01 

1 

2.09*90000  00 

4.4  397  /**(>- oi 

4.600  74  7AU-01 

4.47422770-01 

4.57430360-01 

4.37429090-01 

4 

- J. 341H00O0  0 l 

4. 5? H 7 4600- 01 

4.57316260-01 

4.57429960-01 

*.574 30200-01 

s 

A.UIA/40OU-0I 

4.56*1  1*40-01 

4.57444010-01 

4.5*4  303 10-01 

4.57431710-01 

6 

2.90952000-01 

4. 5 72  76120- CM 

4.57429270-01 

4.574301 70-01 

7 

s.  j4Mioor.-oi 

4.  s  Motion- 01 

4.574  1046U-01 

4.5  7429670-01 

A 

4.21295000-01 

4.57423J7U-01 

4.57430090-01 

9 

4.  nmoou-ot 

4.5  74?AAS0-Q| 

4.57429920- 01 

10 

4 -4052  300U-0 1 

4.57429A6U-01 

11 

4.61  12SOOO-01 

4.57429920-01 

12  4.54105000-01 

I)  4.5*757000-01 
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Table  U.6b  -  Shanks'  Method  Applied  to  the  First  13  Iterations  f  the 
Beal  Matched  Exponent  of  a  Three  Exponential  Approximation  of  f  (t ) 
Using  the  New  Method.  3 


J 

el 

e2 

*3 

eu 

e5 

i.ow.noo 

00 

2. IH5A696U 

00 

7.75445170 

00 

2. >4020960 

00 

2.247141*0 

00 

2.24658140 

00 

7.0655.I0U 

00 

2.29950  7  70 

00 

7.74675500 

00 

2.2^775740 

00 

J.  •  246  7070  J 

OJ 

2.24660690 

00 

00 

2.266925  30 

00 

2.74841990 

00 

2.24730500 

00 

2. 24661600 

00 

2.24660970 

00 

7 .05086500 

00 

2.25496630 

00 

2. 24  ?!  2520 

00 

2. 24720730 

00 

2.24e607?n 

00 

7.457)5401) 

00 

2.25067920 

00 

2.24752020 

00 

2*24606520 

00 

2.24661060 

00 

7.0<»198700 

00 

2 . 24 A  75960 

00 

2.24694550 

00 

2.24703400 

00 

7.58551700 

oo 

2.24793210 

00 

2.2468  7410 

00 

2.24666260 

00 

2*  12953500 

00 

2.24740330 

00 

2.24677060 

On 

2.  ) -.010600 

00 

2.24722420 

00 

2.246674  70 

00 

7. 14854600 

00 

2.24705350 

oo 

2 • i2*00500 

00 

2.24693560 

00 

l?  ?• 100)0000  00 

15  7.50455000  oo 


? .  ?'*r.f,070on  oo 
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V.  DISCUSSION  OF  RESULTS  AND  AREAS  FOR  FURTHER  WORK 

Throughout  this  work  we  have  assumed  the  real  function  f(t)  to 
he:  known  analytically  for  all  time;  piecewise  continuous;  and  of 
bounded  energy,  f^(t)  dt  <  <*>.  Under  these  three  restrictions  we 
reviewed  in  chapter  III  several  ways  of  finding  a  linear  combination 
of  n  exponentials  to  yield  a  least-squares  approximation  of  the  func¬ 
tion  over  the  semi-infinite  interval.  The  results  presented  in  chap¬ 
ter  IV  show  the  new  method  is  the  best  of  these  for  finding  the  matched 
exponents.  This  method  requires  that  the  function  F(s)  and  its  first 
derivative  be  evaluated  only  at  the  n  points  s=  -  in  the  right-half 
of  the  s-plane.  In  the  vicinity  of  these  points  the  function  is  always 
well-behaved,  as  may  be  seen  from  the  Cauchy-Schwartz  inequality, 

OO  00  01 

1  F(— a± )  1 2=  I /  f  (t)exr(s^t)  dt | 2 <_[ /  f2(t)dt][/  lexp^t))2  dt] 
o  o  o 

or 

|F(-s  Ji— 1— .l/V(t)dt]1/2  (5*1} 

/2Re{-s^}  o 

provided  (-s^)  is  in  the  right  half  plane. 

The  restriction  that  the  signal  be  expressed  initially  as  an  analytic 
function  of  time  can  also  be  removed  provided  the  signal  is  expressed 
cn  some  other  basis  such  as  f=  I^k^k  for  which  the  ^(s)  are  known  in 
the  right  half  plane.  The  next  section  gives  an  important  example  in 
which  the  signal  is  represented  initially  on  a  discrete  primal  basis. 

V.l  The  New  Method  Applied  to  Sampled  Data 

Let  p(t)  =  6(t-kT)  denote  an  impulse  train  with  the  impulses 

ft"*® 

spaced  T  second  apart.  The  sampling  of  a  function  can  be  described 


mathematically  by  multiplication  with  p(t),  That  is 

09  m 

f»(t)*p(t)f(t)=  l  f(t)6(t-kT)«  l  f(kT)6(t-kT)  (5.2) 

k=+“  k*-« 

It  is  easily  shown  [33]  that  the  Laplace  transform  of  f*(t)  is 

<X>  00 

F«(s)=  / f»(t)e'8t  dt*  l  f(kt)e"kTa  (5.3) 

o  k=-« 

Consider  the  change  in  variable  z=exp(Ts)  which  maps  the  left  half  plane 
in  the  s  domain  inside  the  unit  circle  in  the  z  domain.  Then  the 
Z-transform  of  the  function  f(t)  (f(t)*0  t  <  0)  is  defined  to  be 


F(z)=  l  f(kT)z  . 
k*o 


(5.U) 


The  approximating  function  at  these  sampled  instants,  is  given  by  the 


rational  Z-transform: 


6  -1  ^  -(n-1) 

a  +a_  z  +...+  a_  z 


- _  /  \  A  T  •  •  *  T  Ok  A 

F  (,)-  44  =  - - — 

avz'  dH)  -1.  .  h  -n 

1+b,  z  +...+  d  z 


°1  t  °2  +  °n 

z  “21  z  ‘Z2  z  •  2n 


(5.5) 


The  poles  of  F  (z)  must  all  be  inside  the  unit  circle  to  ensure  stability 

A 

and  thus  |z^|  >1  k*l,2,...n.  The  error  at  these  sampled  instants 

is  defined  to  be 


e(kT)  ■  f(kT)  -  f  (kT) 

3a 


(5.6) 


E(x)  ■  F(z)  -  F  (z). 

Then  (Ragazzini  p.  179) 

J-  He(kT)]2-  -ij  /  E(z)E(i)  & 
k-o  unit 

circle 


(5.7) 


(5.8) 
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The  necessary  conditions  on  the  2n  parameters  {a^,z^}  to  minimize 
the  functional  J  are 


3J 

-  1  1 

fE(z)  dz 

P 

II 

r  zU-Zj) 

H_.  o 

-  2-i-  i 

rE(z)  dz 

22^J  $ 

r(z-z 

(5.9a) 

(5.9b) 


From  the  Cauchy  integral  formula  and  the  fact  that  E(z)  has  all  its 
poles  inside  the  unit  circle 


E(z.)  -  0 


(5.10) 


E'(*i)  r  o  i=l,2,...,n 

Equations  (5.10)  are  intuitively  correct  Bince  they  are  the  Aigrain- 
Williams  equations  applied  to  sampled  data.  Notice  |z^|  >  1 
corresponds  to  a  point  in  the  right  half  plane  in  the  frequency  domain. 
These  equations  are  solved  iteratively  exactly  as  before  except  one 
uses  F(z)  instead  of  F(s). 

Steiglitz  and  McBride  [3b]  have  also  applied  their  more  complicated 
method  to  sampled  data. 

V.2  Concluding  Remarks 

For  large  n(n>5)  double-precision  arithmetic  is  required  to  get  mean¬ 
ingful  results  using  the  new  method.  This  is  not  unexpected  since  the 
same  difficulty  arises  in  the  simpler  linear  least-squares  approximation 
discussed  in  chapter  II.  Based  on  experience  with  the  tvo  methods,  it 
is  of  the  author's  opinion  that  the  roundoff  errors  in  this  nonlinear 
approximation  will  be  about  the  same  order  of  magnitude  as  those  in 
chapter  II.  The  computational  aspects  of  this  method  deserve  additional 
study,  but  they  will  not  be  pursued  further  in  this  thesis  because  they 
involve  considerations  foreign  to  the  main  thrust  of  this  work. 
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APPENDIX  A 


Construction  of  Orthonoraal  Functions 

Equation  (2.M)  can  also  be  used  in  a  reverse  manner  so  that  if 
JG”^  is  known,  one  can  sometimes  construct  an  orthonormal  basis  by 
simple  inspection  and  avoid  the  Graa-Schmidt  method  altogether*  With 
Gastinel's  result  (derived  without  regard  to  orthogonal  functions) 
and  use  of  (2.1l)  it  is  possible  to  derive  Kautz's  important  result 
for  orthogonalizing  exponentials.  This  second  application  is  demonstrate 
ed  by  the  following  example. 

A  General  Formula  for  Orthonormal  Polynomials  with  Respect  to  a  Con¬ 
stant  Weight  Function 
Define 

s .-1/2 

Xj*  t  1  ,  l/2» 


then 


1  s-1/2  s  -1/2  , 

8ij  ■  /  t  1  t  J  dt  *  y  i„)»l,2 . n.  (A.l) 


As  shewn  previously,  the  inverse  of  this  n  X  n  symmetric  matrix  is 

T,  T.  (A.2) 


«Ij  *  t v^jr  'h 1  j 


where 


TT  V. 
T.  •  1 1  • 
k-i  k  B 
kfta 


Prom  (2.17)  and  (A.2) 


nn 


-1 


2s. 


7?. 


n  n 


(A.  3) 


Also  c  c  .■  g“.,  and  from  (A.2)  and  (A.3) 
nn  nj  n j 

cnj  ”  (2sn)1/2  TT^T  V 


(A.U) 


Since  the  formula  must  hold  for  any  n, 

c  .(2i  )i/2  fr  !£Ll 

cij  t2Bi;  Ts7+s”T  I  I  8  -8 

1  J  k-lJ 

k*4 

Hence  by  (2.3),  (2.U)  and  (A.5) 

1  s-1/2 

♦i(t)  “  l  J  • 

1  4-1  1J 


4  I  i* 


(A.5) 


(A.6) 


However,  this  set  is  orthonormal  on  (0,l).  To  generalize  to  (a,b), 
consider  the  linear  transformation  t»kt'  ♦  d.  When  t«0,t'  -  a  and 
when  t-1,  t'  -  b.  So  k-l/(b-a)  and  d-  -a/(b-a).  Hence,  the  general 


formula  is 


,  n  f(2s  )1; 

♦i(t'Mb-*r  Tv 


'(2a .) 


■ ,“1/2 


4 


(A. 7) 


Note  that  there  is  no  requirement  Sj-l/2  be  an  integer*  Now  let  Pa(t) 
denote  the  Legendre  polynomial  of  degree  n.  It  is  not  hard  to  show 


from  (A.7)  that,  in  this  special  case  for  which  s^-1/2  ■  J-1* 

p  (tw-u"  I  i 

n  4-1  (-2)J”1(n-j) ! ( (4-1)1) 


(A.T  ) 


where 


(A.9) 


+  Note  that  J*  ^(t)  ^(t)  dt  -  /*♦*(*')  ^(t')  dt'/(b-a)-  6^.  Hence, 
the  factor  (b-a)”1^2  appears  in  (A.7). 


+t  The  factor  (-l)n  can  be  dropped  and  the  polynomials  will  still  be  or- 


thonomal.  This  factor  is  added  to  make  (A.8)  agree  with  the  standard 
Legendre  polynomials. 


t  Sxego,  [25]  tec.  11.1.10,  recognised  this  formula  for  orthonormal 
poly  nomi  ala. 


*-r.i ff  &■*■*% 
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Concluding  Remarks 

The  method  described  at  the  beginning  of  chapter  II  shows  a  way 
of  finding  a  closed  form  inverse  of  some  Gram  matrices  that  often  occur 
in  linear  least-squares  problems,  provided  an  analytic  expression  for 
an  appropriate  set  of  orthonormal  functions  can  be  found  in  terms  of 
the  original  basis  elements.  If  an  analytic  expression  cannot  be 
found  for  the  orthonormal  functions,  the  Gram-Schmidt  procedure  can 
always  be  used.  But  then  the  method  loses  seme  of  its  merit,  for  if 
the  basis  elements  are  highly  correlated,  one  may  encounter  the  new 
difficulty  of  computing  the  elements  of  C  accurately. 

Another  distinct  advantage  of  this  method  over  the  "direct''  use 
of  orthonormal  functions  in  least-squares  is  that  it  will  reveal  com¬ 
mon  factors  that  may  be  present  in  each  tem  of  the  inverse.  This  is 
illustrated  by  equation  (B.5)  in  Appendix  B  which  shows  the  common 
factor  [ ( i— 1) I ( J— 1) I ]  of  each  tens  in  the  Inverse  of  G( Laguerre ) . 

It  is  unlikely  that  this  common  factor  would  have  been  observed  if 
linear  combinations  of  the  orthooormal  functions  were  used  to  recon¬ 
struct  the  original  basis.  Finding  such  factors  when  they  exist  can 
obviously  save  time  and  improve  computational  accuracy.  The  results 
for  the  Hilbert  matrix  are  even  better. 

Perhaps  more  important  than  the  direct  application  to  the  least- 
squares  problem  is  the  possibility  of  constructing  orthonormal  beset' 
by  simple  inspection  from  vh*n  •n  •*Plicit  expression 

for  g“*  can  be  found  (as  in  the  case  of  the  Hilbert  matrij.).  The 
construction  of  the  orthooormal  basis  for  fractional  powers  of  t  was 


achieved  by  this  method 
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APPENDIX  B 


A  Least-Squares  Problem  Using  Lsguerre  Polynomials 

The  Laguerre  polynomial  L  (t)  ie  a  polynomial  of  degree  n  in  t 

& 

for  which  M 


/  e-t  L  (t)L  (t)  dt  -  6  . 

DM  DM 

O 


(B.l) 


Leguerre  polynomials  are  orthonormal  with  respect  to  the  weight  function 


e~t  over  (0,  •).  It  ie  known  that 

n  _  ,  /  ^ \k-l 


lo-l(t)  *  klL  O  fell 


(B.2) 


so  that 


•  .<t). 

n— l 


(B.3) 


Suppose  one  desires  to  find  the  such  that  a  continuous  function  f(t) 
is  approximated  in  the  least-square  sense  over  (0,  •)  by 

f  (t)  •  l  a  (tk-Xe“t/2)  -  I  a^x  (t). 

*  k»l  *  k«x  *  * 


(This  may  appear  to  he  sn  odd  choice  of  the  x^,  but  they  are  much  easier 
-t/2 

to  work  vith  than  e  L^( t) ,  Just  as  integrals  involving  single  terms 
of  the  form  expfs^t)  are  easier  to  evaluate  analytically  than  integrals 
involving  orthogonal  functions  formed  from  the  exponentials.)  Then,  as 


in  (2.3)  it  is  immediately  seen  from  (B.2)  that 


(i-u 

lj-r  (j-i)T 


•  o 


i  >  J 
i  <  j. 


(B.M 


Fran  (2.11)  one  finds* 

t  At  first  glance  one  might  think  that  g“X  «■  M  i®  the  Hilbert  ma¬ 

trix.  However,  the  result  held  there  regardless  of  the  ordering  of  the  s^ 
Here  (B.M  holds  only  for  a  particular  ordering  of  the  basis  elements  mid 
the  generalisation  cannot  be  made.  Hsmce,  the  s tarnation  is  necessary. 


>-  V,*.f 


»V"  ■ 


A- 


m 


ft 


where 


end 
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gij*<Vxj>"  /  e-t  ti+J"2  *t  -  (l*j-2)  | 

O 

V  /  rttH^V172  dt  4-1,2,... 

*  o 


(B.5) 


(B.6) 


(B.7) 


Hence,  the  solution  in  closed  fora  is  A  -  g*1?.  (Assigning  thst  (B.7) 
m$v  be  evslunted  in  closed  fora.) 
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APPENDIX  C 


Numerical  Example  for  Mode-1  Iteration 
Consider  the  function 

fu(t)  -  e“ot  [u^t)  -  ujt-l)] 

to  be  approximated  by  one  exponential.  The  Laplace  transform  of  f^(t) 
is 


Fk(s) 


B+O 


and 


(*i ]a 


From  Table  4,  page  61 


where 


and 


V11  gll 
W11  hll 


(a,)4 

X 

i  4 

s 

1 

(bl}j 

yl 

Fu(s) 


ds  1 

2nj  ‘27b ^ 


_  ur  *4'“'  ds  _  U'^lM-l' 

d.  V(bi>.i-1> 


J°  H 

Xi*  "j„,(s+(b1 )  <_1  )(-s^(b1 )  )  2irj  “  ”  2 


BF^Cs) 


rj-r 


l'J-r 


4° 


V®) 


.  wf  -4'"  ds  p4((bl)1-l) 

11 :1- 


j"  PU(8)F4(-»)  ds 

hn‘  J}.  i.*(bl)>1)(-a+lb1,J_1)  2,j  • 

To  evaluate  the  last  integral  by  residues  one  must  be  particularly  care¬ 


ful  because  the  product  F^sjF^-s)  has  an  essential  singularity  at  s*«°. 
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Hen  ce,  one  cannot  make  uae  of  Jordan' a  lemma,  [36]  p.  300,  to  directly 
evaluate  the  integral  by  residues.  Hcwever,  the  integral  may  be  broken 
dcwn  into  the  svsn  of  two  parts,  one  which  vanishes  along  the  infinite 
semi-circular  arc  containing  the  left-half  plane  and  bhe  other  which 
vanishes  along  the  infinite  arc  containing  the  right-half  plane.  That 
is 


J- 

hll"  * 


l+e-2a-e-lB+a) 

(8+a)(s+(b1)^1)(-s+a)(-s+(b1)^1)  2wJ 

RHP 

J*°  -e-(o-s)  d(| 

+  (s+a)(8+(b1)J_1)(-S‘fa)(-sV(b1)j_1)  2*7 

LHP 


dS 


which  simplifies  to 


Finally 


.  l*c"2tt-2e"ta*tbl\l-l* 

11 


1-e 


-2a 


2(bi)J.i(a2-(b1)2.i  )  2a((b1)2.1-a2) 


r  «fu(s)fu(-)  d8 
l"  -  12  "  mL  D.  ,(■)  D,  ,(-■)  2ttJ 


-J-  *>1 


J-l' 


1+e 


-2a 


1+e 


-2a 


2(a2-(bJ,  J  2((bJj_1-a2) 


0. 


rj-r  "  r 

An  iterative  algorithm  very  similar  to  the  one  described  on  page  46 


is  then  used  to  find  (a,),  and  (b, ) , 

X  J  1  J 


81- 


[1]  R. 

[2]  W. 

[3]  T. 

[4]  T. 

[5]  R. 

[6]  L. 

[7]  R. 

[8]  D, 

19]  W 

[10]  0 

Ul]  R 

[12]  P 
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